Load packages

#installation if necessary
#if (!require("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#if (!requireNamespace("devtools", quietly = TRUE)) {
#  install.packages("devtools")
#}

#BiocManager::install("rhdf5")
#BiocManager::install("tximport")
#BiocManager::install("DESeq2")
#BiocManager::install("tximeta")
#install.packages("tidyselect")

#loading in the data
#library(rhdf5) # allows for the handling of hdf5 file formats
#if (!require("tidyverse")) install.packages("tidyverse"); 
#library(tidyverse)
#library(tximport) # package for getting Salmon results into R
#library(datapasta) # allows for copy and pasting
#library(DESeq2)
library(SummarizedExperiment)
library(readr)
#library(tximeta)

#filtering and normalization
#BiocManager::install("vsn")
#library(edgeR)
#library(matrixStats)
#library(cowplot)
#library(hexbin)
#library("vsn")

#multivariate analysis (edgeR)
#if (!require("DT")) install.packages("DT"); library(DT) # for making interactive tables
#if (!require("plotly")) install.packages("plotly"); library(plotly) # for making interactive plots
#if (!require("gt")) install.packages("gt"); library(gt) # A layered 'grammar of tables' - think ggplot, but for tables
#if (!require("stargazer")) install.packages("stargazer"); library(stargazer)
#if (!require("rgl")) install.packages("rgl"); library(rgl)
#if (!require("dendextend")) install.packages("dendextend"); library(dendextend) # A way to customize dendrograms to make more complicated visualizations
#if (!require("ggpattern")) install.packages("ggpattern"); library(ggpattern)
#if (!require("data.table")) install.packages("data.table"); library(data.table)

#making pretty figs and tables
#devtools::install_github("thomasp85/patchwork")
#library(patchwork)
library(gt)
library("ggVennDiagram")
library(ggh4x)
library(patchwork)
library(ggpattern)



#lmerSeq analysis

#install dependencies
#install.packages(pkgs = c('pbapply', 'lmerTest', 'dplyr', 'devtools', 'magrittr', 'gtools', 'nlme', 'numDeriv', 'lavaSearch2'))

#download lmerSeq

#devtools::install_github("stop-pre16/lmerSeq", build_vignettes = T)

#other packages for lmerSeq
#if (!require("DESeq2")) install.packages("DESeq2"); 
#library(DESeq2) 
#if (!require("magrittr")) install.packages("magrittr"); 
#library(magrittr) 
#if (!require("limma")) install.packages("limma"); 
#library(limma)
#if (!require("dendextend")) install.packages("dendenxtend"); 
#library(dendextend) # To create dendrograms
#library(lmerSeq)

#for REBEL analysis
#install dependencies for RESCUE
#scater, scran, scuttle
#BiocManager::install("scater")
#BiocManager::install("scran")
#library(scran) #loads scuttle as part of it

#install Matrix.utils for REBEL
#remotes::install_github("cvarrichio/Matrix.utils") 

#download RESCUE + REBEL - isntructions here: https://github.com/ewynn610/RESCUE and here: https://github.com/ewynn610/REBEL
#devtools::install_github("https://github.com/ewynn610/Rescue",  build_vignettes = T)
#devtools::install_github("https://github.com/ewynn610/REBEL", build_vignettes = TRUE)
library(REBEL)
library(dplyr)
library(ggplot2)
library(plotly)
library(tidyr)

#phylogenies
library(ggtree)
library(ape)
library(treeio)

#general tools
library(stringr)

Clear workspace

rm(list = ls())

Read in data

Gene-level analysis

  • Read in sample metadata
#read in sample metadata
meta <- read.delim("study_design.txt", header = TRUE, as.is = TRUE)

#make the SAMPLE column the rownames
rownames(meta) <- meta$SAMPLE

#add a names column (for lmerSeq later)
meta$names <- meta$SAMPLE

#capture sample labels from the design study file
sampleLabels <- meta$SAMPLE
  • Import variance-stabilized data
#import
vst_expr <- read_tsv("./variance_stabilized_counts_2025_03_12.tsv")

#convert to matrix
expr_mat <- as.matrix(vst_expr[,1:ncol(vst_expr)-1])
rownames(expr_mat) <- vst_expr$geneID
  • Import gene annotations
#get annotations for transcripts by loading file from ncbi ---- WAS EDITED TO MATCH ORs ADDED FROM Mitchell's ANNOTATION----
annotation_file <- ("Ppyr_annotations_ORmodified_2025_03_12.tsv")

#Tx <- read_tsv(annotation_file)
Tx <- read_tsv(annotation_file, col_names = c("target_id", "geneID", "name", "PpyrOR"), col_types = "cccc")

Tx <- as_tibble(Tx)

Analysis: Which genes/ORs are most highly expressed?

  • Calculate average expression across tissue * sex combinations
#calc average gene expression and standard deviation in relevant tissues, gene identity, OR identity
my.pheromone.data.df <- as_tibble(expr_mat) %>% 
  mutate(geneID = row.names(expr_mat)) %>%
  rowwise() %>%
  mutate(antenna.AVG = (SEL534ant + SEL535ant + SEL536ant + SEL543ant + SEL562ant + SEL627ant + SEL672ant)/7,
         antenna.SD = sd(c(SEL534ant, SEL535ant, SEL536ant, SEL543ant, SEL562ant, SEL627ant, SEL672ant)),
         leg.AVG = (SEL534BL + SEL535BL + SEL536BL + SEL543BL + SEL562BL + SEL627BL + SEL672BL)/7, 
         leg.SD = sd(c(SEL534BL, SEL535BL, SEL536BL, SEL543BL, SEL562BL, SEL627BL, SEL672BL)), 
         male.ant.AVG = (SEL534ant + SEL535ant + SEL536ant + SEL543ant)/4,
         male.ant.SD = sd(c(SEL534ant, SEL535ant, SEL536ant, SEL543ant)),
         male.leg.AVG = (SEL534BL + SEL535BL + SEL536BL + SEL543BL)/4,
         male.leg.SD = sd(c(SEL534BL, SEL535BL, SEL536BL, SEL543BL)),
         female.ant.AVG = (SEL562ant + SEL627ant + SEL672ant)/3,
         female.ant.SD = sd(c(SEL562ant, SEL627ant, SEL672ant)),
         female.leg.AVG = (SEL562BL + SEL627BL + SEL672BL)/3,
         female.leg.SD = sd(c(SEL562BL, SEL627BL, SEL672BL)))

#ungroup the rows
my.pheromone.data.df <- my.pheromone.data.df %>%
  ungroup()

#Add annotations (but only distinct geneIDs from the Tx file- has isoforms in it, and want to keep the PpyrOR name associated with the gene ID, even if it only matched 1 isoform)
Tx.with.PpyrORs.per.geneID <- Tx %>% 
  group_by(geneID) %>%
  distinct(PpyrOR, .keep_all = TRUE) %>%
  arrange(PpyrOR)
#checking
#nrow(Tx) #37237
#length(unique(Tx$geneID))#23,919
#nrow(Tx.with.PpyrORs.per.geneID) #23,934

#ungroup
Tx.with.PpyrORs.per.geneID.ungrouped <- Tx.with.PpyrORs.per.geneID %>%
  ungroup()

#get one entry per geneID
Tx.one.entry.per.geneID<- distinct(Tx.with.PpyrORs.per.geneID.ungrouped, geneID, .keep_all = TRUE)
#checking
# nrow(Tx.one.entry.per.geneID) #23,919
# length(unique(Tx.one.entry.per.geneID$PpyrOR)) #103
 
my.pheromone.data.df.expanded <- left_join(my.pheromone.data.df, Tx.one.entry.per.geneID, by = "geneID") 
#nrow(my.pheromone.data.df.expanded) #13,340

#add gene expression ranks by tissue/sex combo
my.pheromone.data.df.expanded$male.ant.RANK <- rank(-my.pheromone.data.df.expanded$male.ant.AVG, ties.method = "min")
my.pheromone.data.df.expanded$male.leg.RANK <- rank(-my.pheromone.data.df.expanded$male.leg.AVG, ties.method = "min")
my.pheromone.data.df.expanded$female.ant.RANK <- rank(-my.pheromone.data.df.expanded$female.ant.AVG, ties.method = "min")
my.pheromone.data.df.expanded$female.leg.RANK <- rank(-my.pheromone.data.df.expanded$female.leg.AVG, ties.method = "min")
#nrow(my.pheromone.data.df.expanded) 13,340

#get a subset of just ORs
my.OR.data.df.expanded <- my.pheromone.data.df.expanded %>%
    filter(!is.na(PpyrOR))

#check to make sure is getting the "correct" subset
#nrow(my.OR.data.df.expanded) #45

### Extract the top 10 genes from each tissue * sex combo

#extract top 10
male.ant.top10 <- my.pheromone.data.df.expanded %>%
  arrange(male.ant.RANK) %>%
  select(name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

#extract top 10
female.ant.top10 <-my.pheromone.data.df.expanded %>%
  arrange(female.ant.RANK) %>%
  select(name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

#extract top 10
male.leg.top10 <- my.pheromone.data.df.expanded %>%
  arrange(male.leg.RANK) %>%
  select(geneID, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

#extract top 10
female.leg.top10 <- my.pheromone.data.df.expanded %>%
  arrange(female.leg.RANK) %>%
  select(geneID, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

Which genes are most highly expressed?

  • What are the top 10 genes in each tissue * sex condition?

Male antennae

#print nice table
male.ant.top10 |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
Gene
Antenna
Leg
M F M F M F M F
PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA 18.97 19.36 1 1 17.99 18.24 8 6
PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA 18.57 18.83 2 2 13.16 13.26 693 576
PREDICTED: Photinus pyralis allergen Tha p 1-like (LOC116176659), mRNA 17.99 18.05 3 3 17.82 17.66 10 12
PREDICTED: Photinus pyralis uncharacterized LOC116179393 (LOC116179393), mRNA 17.86 16.61 4 24 14.10 12.20 309 1,466
PREDICTED: Photinus pyralis probable cytochrome P450 6a14 (LOC116162053), mRNA 17.72 16.96 5 16 8.79 8.48 11,820 13,134
PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA 17.54 17.37 6 8 12.80 12.01 952 1,743
PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA 17.51 17.83 7 4 18.70 18.36 4 4
PREDICTED: Photinus pyralis general odorant-binding protein 19d-like (LOC116167740), mRNA 17.42 17.19 8 11 13.35 13.00 605 734
PREDICTED: Photinus pyralis cytochrome P450 4C1-like (LOC116161249), mRNA 17.20 17.67 9 6 15.66 15.62 75 68
PREDICTED: Photinus pyralis uncharacterized LOC116181535 (LOC116181535), transcript variant X1, mRNA 17.14 17.74 10 5 17.71 17.72 12 10

Female antennae

female.ant.top10 |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
Gene
Antenna
Leg
M F M F M F M F
PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA 18.97 19.36 1 1 17.99 18.24 8 6
PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA 18.57 18.83 2 2 13.16 13.26 693 576
PREDICTED: Photinus pyralis allergen Tha p 1-like (LOC116176659), mRNA 17.99 18.05 3 3 17.82 17.66 10 12
PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA 17.51 17.83 7 4 18.70 18.36 4 4
PREDICTED: Photinus pyralis uncharacterized LOC116181535 (LOC116181535), transcript variant X1, mRNA 17.14 17.74 10 5 17.71 17.72 12 10
PREDICTED: Photinus pyralis cytochrome P450 4C1-like (LOC116161249), mRNA 17.20 17.67 9 6 15.66 15.62 75 68
PREDICTED: Photinus pyralis elongation factor 1-alpha (LOC116166585), transcript variant X1, mRNA 16.69 17.44 17 7 17.16 17.50 19 14
PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA 17.54 17.37 6 8 12.80 12.01 952 1,743
PREDICTED: Photinus pyralis uncharacterized LOC116174483 (LOC116174483), mRNA 16.99 17.22 12 9 16.14 15.84 51 58
PREDICTED: Photinus pyralis general odorant-binding protein 57d (LOC116162318), mRNA 16.54 17.21 22 10 8.84 8.55 11,491 12,964

Male legs

male.leg.top10 |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
geneID Gene
Antenna
Leg
M F M F M F M F
116178108 PREDICTED: Photinus pyralis myosin heavy chain, muscle (LOC116178108), transcript variant X1, mRNA 16.26 16.40 33 33 19.42 19.27 1 2
116180183 PREDICTED: Photinus pyralis actin, muscle (LOC116180183), mRNA 15.90 16.05 46 45 19.33 19.42 2 1
116168454 PREDICTED: Photinus pyralis calcium-transporting ATPase sarcoplasmic/endoplasmic reticulum type (LOC116168454), transcript variant X1, mRNA 16.51 16.37 23 35 19.05 18.55 3 3
116180227 PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA 17.51 17.83 7 4 18.70 18.36 4 4
116173867 PREDICTED: Photinus pyralis paramyosin, long form-like (LOC116173867), mRNA 15.04 14.92 90 110 18.52 18.19 5 7
116174773 PREDICTED: Photinus pyralis titin (LOC116174773), mRNA 14.41 14.53 166 176 18.40 18.25 6 5
116166781 PREDICTED: Photinus pyralis troponin T, skeletal muscle (LOC116166781), transcript variant X1, mRNA 14.93 15.14 99 91 18.01 17.91 7 8
116167561 PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA 18.97 19.36 1 1 17.99 18.24 8 6
116168198 PREDICTED: Photinus pyralis muscle LIM protein Mlp84B-like (LOC116168198), transcript variant X1, mRNA 13.83 13.65 304 399 17.84 17.66 9 11
116176659 PREDICTED: Photinus pyralis allergen Tha p 1-like (LOC116176659), mRNA 17.99 18.05 3 3 17.82 17.66 10 12

Female legs

female.leg.top10 |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
geneID Gene
Antenna
Leg
M F M F M F M F
116180183 PREDICTED: Photinus pyralis actin, muscle (LOC116180183), mRNA 15.90 16.05 46 45 19.33 19.42 2 1
116178108 PREDICTED: Photinus pyralis myosin heavy chain, muscle (LOC116178108), transcript variant X1, mRNA 16.26 16.40 33 33 19.42 19.27 1 2
116168454 PREDICTED: Photinus pyralis calcium-transporting ATPase sarcoplasmic/endoplasmic reticulum type (LOC116168454), transcript variant X1, mRNA 16.51 16.37 23 35 19.05 18.55 3 3
116180227 PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA 17.51 17.83 7 4 18.70 18.36 4 4
116174773 PREDICTED: Photinus pyralis titin (LOC116174773), mRNA 14.41 14.53 166 176 18.40 18.25 6 5
116167561 PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA 18.97 19.36 1 1 17.99 18.24 8 6
116173867 PREDICTED: Photinus pyralis paramyosin, long form-like (LOC116173867), mRNA 15.04 14.92 90 110 18.52 18.19 5 7
116166781 PREDICTED: Photinus pyralis troponin T, skeletal muscle (LOC116166781), transcript variant X1, mRNA 14.93 15.14 99 91 18.01 17.91 7 8
116175784 PREDICTED: Photinus pyralis myosin regulatory light chain 2 (LOC116175784), mRNA 14.71 14.78 118 132 17.78 17.73 11 9
116181535 PREDICTED: Photinus pyralis uncharacterized LOC116181535 (LOC116181535), transcript variant X1, mRNA 17.14 17.74 10 5 17.71 17.72 12 10

How many of the top ten genes are shared?

#display as a Venn Diagram
x = list(male.ant.top10$name, female.ant.top10$name, male.leg.top10$name, female.leg.top10$name)
ggVennDiagram(x, 
              category.names = c("Antenna, M", "Antenna, F", "Leg, M", "Leg, F")) + ggplot2::scale_fill_gradient(low  = "#A0CBE8", high = "#4E79A7")

#sort(intersect(male.ant.top10$name, female.ant.top10$name))
#sort(intersect(male.leg.top10$name, female.leg.top10$name))
#sort(intersect(male.ant.top10$name, male.leg.top10$name))
#sort(intersect(female.ant.top10$name, female.leg.top10$name))
  • Overall, legs look similar to each other and different from antennae - most shared genes are actin, myosin, etc.
  • Two genes are shared among the top 10 lists of all conditions - cytochrome P450 4g15-like (LOC116167561) and ADP,ATP carrier protein (LOC116180227).

Which ORs are most highly expressed?

  • What are the top 10 ORs in each tissue * sex condition?

    • Orco not in top 10 for female legs, though present in all other lists (and first in both male and female antennae).
    • PpyrOR18,23, 92 included in all top 10 lists.
    • PpyrOR6, is highest male antenna-specific (in top 10 list).
    • PpyrOR65 is antenna-specific (in top 10 lists).
    • PpyrOR17PSE in top 10 female legs list (#10) - possibly functional?
### Extract the top 10 ORs from each tissue * sex combo

#extract top 10
male.ant.top10.OR <- my.OR.data.df.expanded %>%
  arrange(male.ant.RANK) %>%
  select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

#extract top 10
female.ant.top10.OR <-my.OR.data.df.expanded %>%
  arrange(female.ant.RANK) %>%
  select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

#extract top 10
male.leg.top10.OR <- my.OR.data.df.expanded %>%
  arrange(male.leg.RANK) %>%
  select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

#extract top 10
female.leg.top10.OR <- my.OR.data.df.expanded %>%
  arrange(female.leg.RANK) %>%
  select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
  slice_head(n = 10)

Male antennae

#print nice table
male.ant.top10.OR |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
PpyrOR Gene
Antenna
Leg
M F M F M F M F
Ppyr/Orco PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA 13.49 13.09 429 670 8.68 8.41 12,528 13,230
PpyrOR6 PREDICTED: Photinus pyralis odorant receptor 94a-like (LOC116160183), transcript variant X1, mRNA 11.84 9.14 1,954 10,155 8.36 8.37 13,291 13,257
PpyrOR18 PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA 9.71 9.79 8,060 7,740 9.13 9.00 10,011 10,647
PpyrOR23 Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA 9.69 10.06 8,149 6,818 9.00 8.98 10,609 10,747
PpyrOR99 PREDICTED: Photinus pyralis odorant receptor 43a-like (LOC116166584), mRNA 9.69 9.14 8,165 10,154 8.49 8.47 13,171 13,153
PpyrOR69 PpyrOR69, NW_022170403 9.43 9.29 9,091 9,549 8.37 8.34 13,285 13,267
PpyrOR92 PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA 9.42 9.84 9,108 7,579 8.71 8.65 12,316 12,593
PpyrOR65 PREDICTED: Photinus pyralis odorant receptor 67c-like (LOC116167462), mRNA 9.36 9.57 9,318 8,499 8.41 8.52 13,259 13,045
PpyrOR100 PREDICTED: Photinus pyralis odorant receptor 43a-like (LOC116166593), mRNA 9.31 8.77 9,544 12,091 8.44 8.50 13,238 13,097
PpyrOR41 PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116175472), mRNA 9.30 9.33 9,578 9,426 8.37 8.38 13,288 13,249

Female antennae

#print nice table
female.ant.top10.OR |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
PpyrOR Gene
Antenna
Leg
M F M F M F M F
Ppyr/Orco PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA 13.49 13.09 429 670 8.68 8.41 12,528 13,230
PpyrOR23 Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA 9.69 10.06 8,149 6,818 9.00 8.98 10,609 10,747
PpyrOR92 PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA 9.42 9.84 9,108 7,579 8.71 8.65 12,316 12,593
PpyrOR18 PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA 9.71 9.79 8,060 7,740 9.13 9.00 10,011 10,647
PpyrOR28 PREDICTED: Photinus pyralis odorant receptor 85c-like (LOC116176545), transcript variant X2, mRNA 9.25 9.67 9,765 8,148 8.87 8.79 11,299 11,781
PpyrOR88 PpyrOR88, NW_022170878 9.19 9.59 10,004 8,445 8.34 8.34 13,297 13,267
PpyrOR65 PREDICTED: Photinus pyralis odorant receptor 67c-like (LOC116167462), mRNA 9.36 9.57 9,318 8,499 8.41 8.52 13,259 13,045
PpyrOR11 PREDICTED: Photinus pyralis uncharacterized LOC116172672 (LOC116172672), mRNA 9.11 9.53 10,271 8,672 8.34 8.34 13,297 13,267
PpyrOR93 PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161161), mRNA 9.05 9.40 10,549 9,158 8.38 8.38 13,278 13,248
PpyrOR70 PpyrOR70, NW_022170403 9.19 9.38 9,983 9,232 8.44 8.34 13,239 13,267

Male legs

#print nice table
male.leg.top10.OR |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
PpyrOR Gene
Antenna
Leg
M F M F M F M F
PpyrOR18 PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA 9.71 9.79 8,060 7,740 9.13 9.00 10,011 10,647
PpyrOR94 PpyrOR94, NW_022170250 9.19 8.95 9,978 11,017 9.01 8.78 10,575 11,844
PpyrOR23 Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA 9.69 10.06 8,149 6,818 9.00 8.98 10,609 10,747
PpyrOR28 PREDICTED: Photinus pyralis odorant receptor 85c-like (LOC116176545), transcript variant X2, mRNA 9.25 9.67 9,765 8,148 8.87 8.79 11,299 11,781
PpyrOR89 PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116180111), transcript variant X1, mRNA 8.84 9.10 11,687 10,292 8.84 9.00 11,483 10,654
PpyrOR9 Similar to XM_031483149.1 PREDICTED: Photinus pyralis odorant receptor Or2-like (LOC116167687), mRNA 8.94 8.99 11,090 10,849 8.82 8.85 11,657 11,448
PpyrOR25 PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116177903), transcript variant X3, mRNA 8.94 8.89 11,091 11,373 8.76 8.81 12,017 11,658
PpyrOR92 PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA 9.42 9.84 9,108 7,579 8.71 8.65 12,316 12,593
Ppyr/Orco PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA 13.49 13.09 429 670 8.68 8.41 12,528 13,230
PpyrOR15 PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116168902), mRNA 9.00 9.05 10,795 10,545 8.67 8.38 12,587 13,251

Female legs

#print nice table
female.leg.top10.OR |>
  gt() |>
  tab_spanner(label = "Antenna", columns = contains("ant")) |>
  tab_spanner(label = "Leg", columns = contains("leg")) |>
  tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
  tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
  cols_label(name = "Gene",
             male.ant.AVG = "M", 
             male.leg.AVG = "M",
             male.ant.RANK = "M",
             male.leg.RANK = "M",
             female.ant.AVG = "F",
             female.leg.AVG = "F",
             female.ant.RANK = "F",
             female.leg.RANK = "F") |>
  fmt_number(columns = ends_with("AVG"),
             decimals = 2) |>
  fmt_number(columns = ends_with("RANK"),
             use_seps = TRUE,
             decimals = 0) |> 
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
  tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
            locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
  tab_style(style = cell_fill(color = "#E5E5E5"),
            locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Antenna", "Leg"))))
Mean Normalized Counts
Rank
Mean Normalized Counts
Rank
PpyrOR Gene
Antenna
Leg
M F M F M F M F
PpyrOR18 PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA 9.71 9.79 8,060 7,740 9.13 9.00 10,011 10,647
PpyrOR89 PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116180111), transcript variant X1, mRNA 8.84 9.10 11,687 10,292 8.84 9.00 11,483 10,654
PpyrOR23 Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA 9.69 10.06 8,149 6,818 9.00 8.98 10,609 10,747
PpyrOR9 Similar to XM_031483149.1 PREDICTED: Photinus pyralis odorant receptor Or2-like (LOC116167687), mRNA 8.94 8.99 11,090 10,849 8.82 8.85 11,657 11,448
PpyrOR25 PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116177903), transcript variant X3, mRNA 8.94 8.89 11,091 11,373 8.76 8.81 12,017 11,658
PpyrOR28 PREDICTED: Photinus pyralis odorant receptor 85c-like (LOC116176545), transcript variant X2, mRNA 9.25 9.67 9,765 8,148 8.87 8.79 11,299 11,781
PpyrOR94 PpyrOR94, NW_022170250 9.19 8.95 9,978 11,017 9.01 8.78 10,575 11,844
PpyrOR92 PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA 9.42 9.84 9,108 7,579 8.71 8.65 12,316 12,593
PpyrOR19 PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178510), transcript variant X2, mRNA 8.72 8.95 12,455 11,018 8.60 8.57 12,866 12,912
PpyrOR17PSE PpyrOR17PSE, NW_022170819 8.55 8.66 13,093 12,641 8.53 8.53 13,075 13,020

How many of the top ten ORs are shared?

#display as a Venn Diagram
top.10.OR.lists <- list(male.ant.top10.OR$PpyrOR, female.ant.top10.OR$PpyrOR, male.leg.top10.OR$PpyrOR, female.leg.top10.OR$PpyrOR)
ggVennDiagram(top.10.OR.lists, 
              category.names = c("Antenna, M", "Antenna, F", "Leg, M", "Leg, F")) + ggplot2::scale_fill_gradient(low  = "#A0CBE8", high = "#4E79A7")

Overall, legs and antennae look different from each other, though sex difference is greater in antennae. If the ORs are shared between M + F antennae, then they are also usually shared with legs.


Which ORs are NOT expressed?

ORs.only <- Tx[!is.na(Tx$PpyrOR),]

#nrow(ORs.only) #103
ORs.not.expressed <- sort(setdiff(ORs.only$PpyrOR, my.OR.data.df.expanded$PpyrOR))

n_PSE <- length(grep("PSE", Tx$PpyrOR)) #15 PSEs
n_PSE.not.expressed <- length(grep("PSE", ORs.not.expressed)) #12 PSEs
  • 58 (56%) ORs are “not expressed” - they didn’t pass filtering in the expression dataset. They may be expressed in other life stages.
  • 12 of the 15 (80%) pseudogenized ORs (PSE) are in the “not expressed” list.
  • 3 (20%) OR PSEs did pass filtering and are expressed in adults, either in antennae or legs.

Analysis: Differential expression

  • Fit the model
#set factors for analysis
meta$TISSUE = factor(meta$TISSUE, levels = c("leg", "antenna"))
meta$SEX = factor(meta$SEX, levels = c("F", "M"))

#convert to Summarized Experiment format
rse = SummarizedExperiment(assays = SimpleList(normcounts = as.matrix(expr_mat)), colData = meta)

#fit the model
res_reb = rebelFit(object = rse, fixedEffects = ~ TISSUE * SEX, #takes awhile
                   subjectVariable = "INDIVIDUAL",
                   pseudoBulk = T, 
                   parallel = F,
                   outputFits = TRUE)  
  • Calculate contrasts
#calculate DE for all contrasts
#male tissues: antenna vs leg
reb.mant_v_mleg <- rebelTest(res_reb, contrast = c(0, 1, 0, 1)) #male antennae vs leg
colnames(reb.mant_v_mleg) <- c("Estimate_mant_v_mleg", "Std.Error_mant_v_mleg", "t.value_mant_v_mleg", "df_mant_v_mleg", "p_val_raw_mant_v_mleg", "p_val_adj_mant_v_mleg")
reb.mant_v_mleg$geneID <- row.names(reb.mant_v_mleg) #add row names

#female tissues: antenna vs leg
reb.fant_v_fleg = rebelTest(res_reb, contrast = c(0, 1, 0, 0)) #female antennae vs leg - same as lmerseq?
colnames(reb.fant_v_fleg) <- c("Estimate_fant_v_fleg", "Std.Error_fant_v_fleg", "t.value_fant_v_fleg", "df_fant_v_fleg", "p_val_raw_fant_v_fleg", "p_val_adj_fant_v_fleg")
reb.fant_v_fleg$geneID <- row.names(reb.fant_v_fleg) #add row names

#antenna: male vs female
reb.mant_v_fant = rebelTest(res_reb, contrast = c(0, 0, 1, 1)) #male antennae vs female antenna - same as lmerseq?
colnames(reb.mant_v_fant) <- c("Estimate_mant_v_fant", "Std.Error_mant_v_fant", "t.value_mant_v_fant", "df_mant_v_fant", "p_val_raw_mant_v_fant", "p_val_adj_mant_v_fant")
reb.mant_v_fant$geneID <- row.names(reb.mant_v_fant) #add row names

#leg: male vs female
reb.mleg_v_fleg = rebelTest(res_reb, contrast = c(0, 0, 1, 0)) #male leg vs female leg - same as lmerseq?
colnames(reb.mleg_v_fleg) <- c("Estimate_mleg_v_fleg", "Std.Error_mleg_v_fleg", "t.value_mleg_v_fleg", "df_mleg_v_fleg", "p_val_raw_mleg_v_fleg", "p_val_adj_mleg_v_fleg")
reb.mleg_v_fleg$geneID <- row.names(reb.mleg_v_fleg) #add row names

#get into 1 big data table
#are all gene IDs in the same order? if so - can just cbind these together? #or just do join_by?
reb.tissue.comp <- left_join(reb.mant_v_mleg, reb.fant_v_fleg, by = "geneID")
reb.sex.comp <- left_join(reb.mant_v_fant, reb.mleg_v_fleg, by = "geneID")
reb.comp <- left_join(reb.tissue.comp, reb.sex.comp, by = "geneID")

#combine with my.pheromone.data.df.expanded (has annotations)
my.DE.data.df.expanded <- left_join(my.pheromone.data.df.expanded, reb.comp, by = "geneID")

#get ORs just on own
my.DE.OR.data.df.expanded <- my.DE.data.df.expanded %>%
  filter(!is.na(PpyrOR))

#nrow(my.DE.OR.data.df.expanded) #45 - yes!
#export table (in nice order) as Supp table (big table)
#write.table(my.DE.data.df.expanded, "2025_03_12_DE_table_expanded.txt", row.names = FALSE, sep ='\t')
  • Collate and organize results
my.DE.data.df.expanded <- my.DE.data.df.expanded %>%
  mutate(OR = case_when(!is.na(PpyrOR) ~ "OR",
                       TRUE ~"Other"),
         DE_mant_v_mleg = ifelse(p_val_adj_mant_v_mleg < 0.05, "Significant", "Not significant"),
         DE_fant_v_fleg = ifelse(p_val_adj_fant_v_fleg < 0.05, "Significant", "Not significant"),
         DE_mant_v_fant = ifelse(p_val_adj_mant_v_fant < 0.05, "Significant", "Not significant"),
         DE_mleg_v_fleg = ifelse(p_val_adj_mleg_v_fleg < 0.05, "Significant", "Not significant"),
         Direction_mant_v_mleg = case_when(
 #          Estimate_mant_v_mleg > 0.6 & p_val_adj_mant_v_mleg < 0.05 ~ "Upregulated",
 #          Estimate_mant_v_mleg < -0.6 & p_val_adj_mant_v_mleg < 0.05 ~ "Downregulated",
           Estimate_mant_v_mleg > 0 & p_val_adj_mant_v_mleg < 0.05 ~ "Upregulated",
           Estimate_mant_v_mleg < 0 & p_val_adj_mant_v_mleg < 0.05 ~ "Downregulated",
           TRUE ~ "Not significant"),
         Direction_fant_v_fleg = case_when(
#           Estimate_fant_v_fleg > 0.6 & p_val_adj_fant_v_fleg < 0.05 ~ "Upregulated",
#           Estimate_fant_v_fleg < -0.6 & p_val_adj_fant_v_fleg < 0.05 ~ "Downregulated",
           Estimate_fant_v_fleg > 0 & p_val_adj_fant_v_fleg < 0.05 ~ "Upregulated",
           Estimate_fant_v_fleg < 0 & p_val_adj_fant_v_fleg < 0.05 ~ "Downregulated",
           TRUE ~ "Not significant"),
        Direction_mant_v_fant = case_when(
#           Estimate_mant_v_fant > 0.6 & p_val_adj_mant_v_fant < 0.05 ~ "Upregulated",
#           Estimate_mant_v_fant < -0.6 & p_val_adj_mant_v_fant < 0.05 ~ "Downregulated",
           Estimate_mant_v_fant > 0 & p_val_adj_mant_v_fant < 0.05 ~ "Upregulated",
           Estimate_mant_v_fant < 0 & p_val_adj_mant_v_fant < 0.05 ~ "Downregulated",
           TRUE ~ "Not significant"),
        Direction_mleg_v_fleg = case_when(
#           Estimate_mleg_v_fleg > 0.6 & p_val_adj_mleg_v_fleg < 0.05 ~ "Upregulated",
#           Estimate_mleg_v_fleg < -0.6 & p_val_adj_mleg_v_fleg < 0.05 ~ "Downregulated",
            Estimate_mleg_v_fleg > 0 & p_val_adj_mleg_v_fleg < 0.05 ~ "Upregulated",
           Estimate_mleg_v_fleg < 0 & p_val_adj_mleg_v_fleg < 0.05 ~ "Downregulated",
           TRUE ~ "Not significant"))

my.DE.data.df.expanded$Direction_mant_v_mleg <- factor(my.DE.data.df.expanded$Direction_mant_v_mleg, levels = c("Upregulated", "Not significant", "Downregulated"))
my.DE.data.df.expanded$Direction_fant_v_fleg <- factor(my.DE.data.df.expanded$Direction_fant_v_fleg, levels = c("Upregulated", "Not significant", "Downregulated"))
my.DE.data.df.expanded$Direction_mant_v_fant <- factor(my.DE.data.df.expanded$Direction_mant_v_fant, levels = c("Upregulated", "Not significant", "Downregulated"))
my.DE.data.df.expanded$Direction_mleg_v_fleg <- factor(my.DE.data.df.expanded$Direction_mleg_v_fleg, levels = c("Upregulated", "Not significant", "Downregulated"))

#the OR table
my.DE.OR.data.df.expanded <- my.DE.data.df.expanded %>%
  filter(!is.na(PpyrOR))

#Orco table
my.DE.ORco.data.df.expanded <- my.DE.data.df.expanded %>%
  filter(PpyrOR == "Ppyr/Orco")

Contrasts

Male: antenna v leg

How many genes and ORs are DE in male antenna vs leg?

m_tissue_DE <- my.DE.data.df.expanded %>%
  filter(DE_mant_v_mleg == "Significant")

mant_v_mleg_DE_ORs <- my.DE.OR.data.df.expanded %>%
  filter(DE_mant_v_mleg == "Significant") %>%
  arrange(desc(Estimate_mant_v_mleg), p_val_adj_mant_v_mleg) %>%
  select(PpyrOR, Direction_mant_v_mleg, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
  • 3678 genes are significantly differentially expressed between male antenna and legs
  • 32 ORs, out of 45 that passed the filter, (71.11)% are DE between male antenna and leg.
Which genes and ORs are they?
#for ggplot
my.mant_v_mleg.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), color = Direction_mant_v_mleg)) +
#  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) + 
  labs(color = "Expression",
       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated"))

#add circles, label Orco, OR6, male-specific ORs, and female-specific ORs
my.mant_v_mleg.plot.volcano.labels <- my.mant_v_mleg.plot.volcano +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_mleg, y=-log10(p_val_adj_mant_v_mleg), fill = Direction_mant_v_mleg), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#F28E2B", "#989ca3", "#4E79A7")) +
  guides(fill = "none") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
  #DE in male antennae
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
  #male-specific DE
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), label = "OR100", color = "black", size = 3, hjust = -1.2, vjust = -1) +
    geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.15, yend = -log10(p_val_adj_mant_v_mleg) + .18), color = "grey28") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), label = "OR67", color = "black", size = 3, hjust = -1.5, vjust = -0.9) +
    geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.7, yend = -log10(p_val_adj_mant_v_mleg) + 0.18), color = "grey28") +
  #female-specific # F: 23, 28, 81, 15, 64, 59PSE

  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), label = "OR59PSE", color = "black", size =3, hjust = -0.88, vjust = 1.1) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.8, yend = -log10(p_val_adj_mant_v_mleg) + -0.1), color = "grey28") +
  
    geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), label = "OR23", color = "black", size = 3, hjust = -1.35, vjust = 1.37) +
    geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.4, yend = -log10(p_val_adj_mant_v_mleg) + -0.14), color = "grey28") +

  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), label = "OR28", color = "black", size =3, hjust = -1.52, vjust = 1.6) +
      geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.73, yend = -log10(p_val_adj_mant_v_mleg) + -0.17), color = "grey28") +
  
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), label = "OR15", color = "black", size = 3, hjust = -1.56, vjust = 1.82) +
    geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.8, yend = -log10(p_val_adj_mant_v_mleg) + -0.2), color = "grey28") +
  
    geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), label = "OR81", color = "black", size = 3, hjust = -1.65, vjust = 2.2) +
    geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.95, yend = -log10(p_val_adj_mant_v_mleg) + -0.25), color = "grey28") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), label = "OR64", color = "black", size = 3, hjust = -1.64, vjust = 3) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.89, yend = -log10(p_val_adj_mant_v_mleg) + -0.35), color = "grey28") +  
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

my.mant_v_mleg.plot.volcano.labels

#for plotly
my.mant_v_mleg.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), color = Direction_mant_v_mleg, text = name)) +
  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) + 
  labs(color = "Male antenna v leg",
#       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
        x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?

my.mant_v_mleg.plot.circles <- my.mant_v_mleg.plot +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_mleg, y=-log10(p_val_adj_mant_v_mleg), fill = Direction_mant_v_mleg, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#F28E2B", "#989ca3", "#4E79A7")) +
  guides(fill = "none")

ggplotly(my.mant_v_mleg.plot.circles)

Circle with bold outline = ORs that passed filtering


Top 25 DE gene list
mant_v_mleg_DE_genes <- my.DE.data.df.expanded %>%
  filter(p_val_adj_mant_v_mleg < 0.05) %>%
  arrange(desc(Estimate_mant_v_mleg), p_val_adj_mant_v_mleg) %>%
  select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)

#nrow(mant_v_mleg_DE_genes) #3678 genes are DE

mant_v_mleg_DE_genes.top25 <- mant_v_mleg_DE_genes %>%
  slice_head(n = 25)

mant_v_mleg_DE_genes.top25 |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(name = "Gene",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
Gene
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p
PREDICTED: Photinus pyralis probable cytochrome P450 6a14 (LOC116162053), mRNA 8.93 0.00 8.48 0.00 0.77 0.23 0.31 0.46
PREDICTED: Photinus pyralis general odorant-binding protein 57d (LOC116162318), mRNA 7.69 0.00 8.66 0.00 −0.68 0.13 0.29 0.37
PREDICTED: Photinus pyralis protein takeout-like (LOC116181411), mRNA 7.64 0.00 8.36 0.00 −0.58 0.29 0.14 0.72
PREDICTED: Photinus pyralis ejaculatory bulb-specific protein 3-like (LOC116176666), mRNA 7.55 0.00 7.22 0.00 0.73 0.26 0.40 0.39
PREDICTED: Photinus pyralis lipase 3-like (LOC116180894), transcript variant X1, mRNA 7.48 0.00 7.45 0.00 0.12 0.82 0.09 0.73
PREDICTED: Photinus pyralis uncharacterized LOC116165182 (LOC116165182), mRNA 7.47 0.00 7.56 0.00 0.07 0.96 0.15 0.72
PREDICTED: Photinus pyralis uncharacterized LOC116172495 (LOC116172495), mRNA 7.37 0.00 7.14 0.00 −0.36 0.51 −0.59 0.16
PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116170590), mRNA 6.61 0.00 6.73 0.00 0.06 0.97 0.18 0.73
PREDICTED: Photinus pyralis uncharacterized LOC116179743 (LOC116179743), transcript variant X1, mRNA 6.28 0.00 4.95 0.00 1.74 0.11 0.41 0.54
PREDICTED: Photinus pyralis general odorant-binding protein 72-like (LOC116170373), mRNA 6.00 0.00 6.34 0.00 −0.42 0.63 −0.09 0.88
PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA 5.41 0.00 5.57 0.00 −0.26 0.72 −0.10 0.83
PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161708), mRNA 5.35 0.00 6.00 0.00 −0.11 0.93 0.53 0.28
PREDICTED: Photinus pyralis farnesol dehydrogenase-like (LOC116162028), mRNA 5.32 0.00 5.17 0.00 0.25 0.75 0.10 0.83
PREDICTED: Photinus pyralis probable cytochrome P450 6a13 (LOC116161618), mRNA 5.00 0.00 4.91 0.00 0.12 0.96 0.02 0.97
PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161490), transcript variant X1, mRNA 4.93 0.00 4.38 0.00 1.12 0.50 0.57 0.55
PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA 4.82 0.00 4.68 0.00 0.40 0.51 0.27 0.46
PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116164250), mRNA 4.79 0.00 4.67 0.00 0.51 0.48 0.39 0.40
PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA 4.75 0.00 5.36 0.00 0.17 0.94 0.79 0.32
PREDICTED: Photinus pyralis uncharacterized LOC116180967 (LOC116180967), mRNA 4.30 0.00 4.35 0.00 0.22 0.90 0.27 0.69
PREDICTED: Photinus pyralis odorant-binding protein 59a (LOC116179853), mRNA 4.23 0.00 5.41 0.00 0.04 0.98 1.22 0.11
PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161159), transcript variant X1, mRNA 4.20 0.00 3.74 0.00 1.38 0.28 0.92 0.33
PREDICTED: Photinus pyralis uncharacterized LOC116178615 (LOC116178615), mRNA 4.20 0.00 4.67 0.00 0.21 0.86 0.69 0.28
PREDICTED: Photinus pyralis glutathione S-transferase 1-like (LOC116167738), mRNA 4.08 0.00 4.17 0.00 0.33 0.85 0.42 0.59
PREDICTED: Photinus pyralis general odorant-binding protein 19d-like (LOC116167740), mRNA 4.08 0.00 4.19 0.00 0.23 0.83 0.35 0.49
PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116176983), transcript variant X1, mRNA 3.97 0.00 4.20 0.00 −0.04 0.98 0.19 0.73

DE OR list
mant_v_mleg_DE_ORs %>%
  select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg) |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(PpyrOR = "OR",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
OR
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p
Ppyr/Orco 4.82 0.00 4.68 0.00 0.40 0.51 0.27 0.46
PpyrOR6 3.48 0.00 0.77 0.03 2.71 0.00 −0.01 0.99
PpyrOR99 1.20 0.00 0.67 0.01 0.55 0.13 0.02 0.95
PpyrOR69 1.05 0.00 0.95 0.01 0.14 0.81 0.03 0.92
PpyrOR65 0.95 0.01 1.05 0.01 −0.21 0.69 −0.11 0.73
PpyrOR41 0.92 0.00 0.95 0.00 −0.03 0.96 −0.01 0.98
PpyrOR100 0.87 0.00 0.27 0.07 0.54 0.09 −0.06 0.76
PpyrOR88 0.84 0.01 1.25 0.00 −0.40 0.28 0.00 1.00
PpyrOR11 0.77 0.01 1.19 0.01 −0.42 0.30 0.00 1.00
PpyrOR70 0.75 0.01 1.04 0.01 −0.19 0.76 0.10 0.79
PpyrOR38 0.73 0.00 0.67 0.00 0.06 0.90 0.00 1.00
PpyrOR62 0.72 0.01 0.69 0.02 0.05 0.96 0.02 0.95
PpyrOR92 0.71 0.02 1.19 0.01 −0.42 0.41 0.06 0.87
PpyrOR63 0.70 0.01 0.85 0.01 −0.18 0.68 −0.03 0.92
PpyrOR93 0.67 0.00 1.02 0.00 −0.35 0.19 0.00 1.00
PpyrOR75 0.63 0.01 0.61 0.01 0.03 0.97 0.01 0.97
PpyrOR24 0.60 0.02 0.99 0.01 −0.42 0.32 −0.03 0.94
PpyrOR71 0.59 0.00 0.72 0.01 −0.07 0.90 0.05 0.81
PpyrOR68 0.59 0.01 0.48 0.02 0.20 0.58 0.09 0.69
PpyrOR56PSE 0.58 0.04 0.69 0.04 −0.09 0.94 0.02 0.97
PpyrOR18 0.58 0.01 0.80 0.01 −0.08 0.92 0.13 0.66
PpyrOR16 0.58 0.01 0.82 0.01 −0.22 0.56 0.03 0.92
PpyrOR98 0.57 0.02 0.85 0.01 −0.26 0.68 0.02 0.97
PpyrOR12 0.56 0.00 0.64 0.01 −0.10 0.78 −0.03 0.90
PpyrOR76 0.55 0.01 0.77 0.00 −0.30 0.42 −0.08 0.75
PpyrOR53 0.52 0.01 1.00 0.00 −0.51 0.11 −0.03 0.91
PpyrOR58 0.49 0.02 0.59 0.02 −0.07 0.93 0.03 0.91
PpyrOR77 0.44 0.01 0.66 0.01 −0.15 0.65 0.06 0.77
PpyrOR61 0.42 0.01 0.38 0.03 0.04 0.96 0.00 1.00
PpyrOR67 0.38 0.04 0.42 0.05 −0.06 0.96 −0.02 0.96
PpyrOR66 0.35 0.02 0.51 0.01 −0.24 0.36 −0.08 0.64
PpyrOR22 0.30 0.02 0.61 0.01 −0.29 0.28 0.02 0.93

Female: antenna v leg

How many genes and ORs are DE in female antenna vs leg?

f_tissue_DE <- my.DE.data.df.expanded %>%
  filter(DE_fant_v_fleg == "Significant")

fant_v_fleg_DE_ORs <- my.DE.OR.data.df.expanded %>%
  filter(DE_fant_v_fleg == "Significant") %>%
  arrange(desc(Estimate_fant_v_fleg), p_val_adj_fant_v_fleg) %>%
  select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
  • 3027 genes are significantly differentially expressed between female antenna and legs
  • 36 ORs, out of 45 that passed the filter, (80)% are DE between female antenna and leg.
Which genes and ORs are they?
#for ggplot
my.fant_v_fleg.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), color = Direction_fant_v_fleg)) +
#  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) + 
  labs(color = "Female antenna v leg",
       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated"))

#add circles, label Orco, OR6, OR67 and 100, OR 15,23,28,59PSE,64,81
my.fant_v_fleg.plot.volcano.labels <- my.fant_v_fleg.plot.volcano +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_fant_v_fleg, y=-log10(p_val_adj_fant_v_fleg), fill = Direction_fant_v_fleg), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#F28E2B", "#989ca3", "#4E79A7")) +
  guides(fill = "none") +
  theme(legend.position = "none") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.05, vjust = -0.45) +
  
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), label = "OR81", color = "black", size = 3, hjust = -2.22, vjust = -1.8) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4, yend = -log10(p_val_adj_fant_v_fleg) + 0.3), color = "grey28") +
  
   geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), label = "OR28", color = "black", size = 3, hjust = -2.17, vjust = -3.6) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 3.9, yend = -log10(p_val_adj_fant_v_fleg) + 0.55), color = "grey28") +
  
   geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), label = "OR59PSE", color = "black", size = 3, hjust = -1.35, vjust = -2.9) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.26, yend = -log10(p_val_adj_fant_v_fleg) + 0.45), color = "grey28") +
   
   geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -2.85, vjust = -2.4) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.05, yend = -log10(p_val_adj_fant_v_fleg) +  0.39), color = "grey28") +
  
     geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), label = "OR64", color = "black", size = 3, hjust = -2.37, vjust = -1.8) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.25, yend = -log10(p_val_adj_fant_v_fleg) + 0.315), color = "grey28") +
  
    geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), label = "OR23", color = "black", size = 3, hjust = -2.07, vjust = -0.6) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 3.7, yend = -log10(p_val_adj_fant_v_fleg) + 0.132), color = "grey28") +
  
    geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), label = "OR15", color = "black", size = 3, hjust = -2.3, vjust = 0.4) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.15, yend = -log10(p_val_adj_fant_v_fleg) + 0.07), color = "grey28") +
  
   geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), label = "OR67", color = "black", size = 3, hjust = -2.42, vjust = 1.4) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.36, yend = -log10(p_val_adj_fant_v_fleg) + -0.12), color = "grey28") +
  
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), label = "OR100", color = "black", size = 3, hjust = -2.07, vjust = 1.78) +
    geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.5, yend = -log10(p_val_adj_fant_v_fleg) + -0.19), color = "grey28") +
  
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))


my.fant_v_fleg.plot.volcano.labels

#for plotly
my.fant_v_fleg.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), color = Direction_fant_v_fleg, text = name)) +
  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) + 
  labs(color = "Female antenna v leg",
#       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
        x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?

my.fant_v_fleg.plot.circles <- my.fant_v_fleg.plot +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_fant_v_fleg, y=-log10(p_val_adj_fant_v_fleg), fill = Direction_fant_v_fleg, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#F28E2B", "#989ca3", "#4E79A7")) +
  guides(fill = "none")

ggplotly(my.fant_v_fleg.plot.circles)

Circle with bold outline = ORs that passed filtering


Top 25 DE gene list
fant_v_fleg_DE_genes <- my.DE.data.df.expanded %>%
  filter(p_val_adj_fant_v_fleg < 0.05) %>%
  arrange(desc(Estimate_fant_v_fleg), p_val_adj_fant_v_fleg) %>%
  select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)

fant_v_fleg_DE_genes.top25 <- fant_v_fleg_DE_genes %>%
  slice_head(n = 25)

fant_v_fleg_DE_genes.top25 |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(name = "Gene",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
Gene
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p
PREDICTED: Photinus pyralis general odorant-binding protein 57d (LOC116162318), mRNA 7.69 0.00 8.66 0.00 −0.68 0.13 0.29 0.37
PREDICTED: Photinus pyralis probable cytochrome P450 6a14 (LOC116162053), mRNA 8.93 0.00 8.48 0.00 0.77 0.23 0.31 0.46
PREDICTED: Photinus pyralis protein takeout-like (LOC116181411), mRNA 7.64 0.00 8.36 0.00 −0.58 0.29 0.14 0.72
PREDICTED: Photinus pyralis uncharacterized LOC116165182 (LOC116165182), mRNA 7.47 0.00 7.56 0.00 0.07 0.96 0.15 0.72
PREDICTED: Photinus pyralis lipase 3-like (LOC116180894), transcript variant X1, mRNA 7.48 0.00 7.45 0.00 0.12 0.82 0.09 0.73
PREDICTED: Photinus pyralis ejaculatory bulb-specific protein 3-like (LOC116176666), mRNA 7.55 0.00 7.22 0.00 0.73 0.26 0.40 0.39
PREDICTED: Photinus pyralis uncharacterized LOC116172495 (LOC116172495), mRNA 7.37 0.00 7.14 0.00 −0.36 0.51 −0.59 0.16
PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116170590), mRNA 6.61 0.00 6.73 0.00 0.06 0.97 0.18 0.73
PREDICTED: Photinus pyralis general odorant-binding protein 72-like (LOC116170373), mRNA 6.00 0.00 6.34 0.00 −0.42 0.63 −0.09 0.88
PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161708), mRNA 5.35 0.00 6.00 0.00 −0.11 0.93 0.53 0.28
PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA 5.41 0.00 5.57 0.00 −0.26 0.72 −0.10 0.83
PREDICTED: Photinus pyralis odorant-binding protein 59a (LOC116179853), mRNA 4.23 0.00 5.41 0.00 0.04 0.98 1.22 0.11
PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA 4.75 0.00 5.36 0.00 0.17 0.94 0.79 0.32
PREDICTED: Photinus pyralis farnesol dehydrogenase-like (LOC116162028), mRNA 5.32 0.00 5.17 0.00 0.25 0.75 0.10 0.83
PREDICTED: Photinus pyralis uncharacterized LOC116179743 (LOC116179743), transcript variant X1, mRNA 6.28 0.00 4.95 0.00 1.74 0.11 0.41 0.54
PREDICTED: Photinus pyralis probable cytochrome P450 6a13 (LOC116161618), mRNA 5.00 0.00 4.91 0.00 0.12 0.96 0.02 0.97
PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA 4.82 0.00 4.68 0.00 0.40 0.51 0.27 0.46
PREDICTED: Photinus pyralis uncharacterized LOC116178615 (LOC116178615), mRNA 4.20 0.00 4.67 0.00 0.21 0.86 0.69 0.28
PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116164250), mRNA 4.79 0.00 4.67 0.00 0.51 0.48 0.39 0.40
PREDICTED: Photinus pyralis uncharacterized LOC116172452 (LOC116172452), mRNA 3.69 0.00 4.43 0.00 −0.54 0.62 0.19 0.77
PREDICTED: Photinus pyralis uncharacterized LOC116179393 (LOC116179393), mRNA 3.76 0.01 4.41 0.01 1.25 0.47 1.90 0.17
PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161490), transcript variant X1, mRNA 4.93 0.00 4.38 0.00 1.12 0.50 0.57 0.55
PREDICTED: Photinus pyralis uncharacterized LOC116180967 (LOC116180967), mRNA 4.30 0.00 4.35 0.00 0.22 0.90 0.27 0.69
PREDICTED: Photinus pyralis sensory neuron membrane protein 1-like (LOC116168667), mRNA 3.94 0.00 4.26 0.00 −0.14 0.92 0.18 0.71
PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116176983), transcript variant X1, mRNA 3.97 0.00 4.20 0.00 −0.04 0.98 0.19 0.73

DE OR list
fant_v_fleg_DE_ORs |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(PpyrOR = "OR",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
OR
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p
Ppyr/Orco 4.82 0.00 4.68 0.00 0.40 0.51 0.27 0.46
PpyrOR88 0.84 0.01 1.25 0.00 −0.40 0.28 0.00 1.00
PpyrOR92 0.71 0.02 1.19 0.01 −0.42 0.41 0.06 0.87
PpyrOR11 0.77 0.01 1.19 0.01 −0.42 0.30 0.00 1.00
PpyrOR23 0.69 0.07 1.08 0.04 −0.37 0.81 0.02 0.98
PpyrOR65 0.95 0.01 1.05 0.01 −0.21 0.69 −0.11 0.73
PpyrOR70 0.75 0.01 1.04 0.01 −0.19 0.76 0.10 0.79
PpyrOR93 0.67 0.00 1.02 0.00 −0.35 0.19 0.00 1.00
PpyrOR53 0.52 0.01 1.00 0.00 −0.51 0.11 −0.03 0.91
PpyrOR24 0.60 0.02 0.99 0.01 −0.42 0.32 −0.03 0.94
PpyrOR69 1.05 0.00 0.95 0.01 0.14 0.81 0.03 0.92
PpyrOR41 0.92 0.00 0.95 0.00 −0.03 0.96 −0.01 0.98
PpyrOR28 0.38 0.10 0.88 0.02 −0.42 0.45 0.08 0.85
PpyrOR63 0.70 0.01 0.85 0.01 −0.18 0.68 −0.03 0.92
PpyrOR98 0.57 0.02 0.85 0.01 −0.26 0.68 0.02 0.97
PpyrOR16 0.58 0.01 0.82 0.01 −0.22 0.56 0.03 0.92
PpyrOR18 0.58 0.01 0.80 0.01 −0.08 0.92 0.13 0.66
PpyrOR76 0.55 0.01 0.77 0.00 −0.30 0.42 −0.08 0.75
PpyrOR6 3.48 0.00 0.77 0.03 2.71 0.00 −0.01 0.99
PpyrOR81 0.17 0.18 0.77 0.01 −0.63 0.09 −0.03 0.90
PpyrOR71 0.59 0.00 0.72 0.01 −0.07 0.90 0.05 0.81
PpyrOR56PSE 0.58 0.04 0.69 0.04 −0.09 0.94 0.02 0.97
PpyrOR62 0.72 0.01 0.69 0.02 0.05 0.96 0.02 0.95
PpyrOR38 0.73 0.00 0.67 0.00 0.06 0.90 0.00 1.00
PpyrOR15 0.33 0.14 0.67 0.04 −0.05 0.98 0.29 0.63
PpyrOR99 1.20 0.00 0.67 0.01 0.55 0.13 0.02 0.95
PpyrOR77 0.44 0.01 0.66 0.01 −0.15 0.65 0.06 0.77
PpyrOR12 0.56 0.00 0.64 0.01 −0.10 0.78 −0.03 0.90
PpyrOR22 0.30 0.02 0.61 0.01 −0.29 0.28 0.02 0.93
PpyrOR75 0.63 0.01 0.61 0.01 0.03 0.97 0.01 0.97
PpyrOR58 0.49 0.02 0.59 0.02 −0.07 0.93 0.03 0.91
PpyrOR64 0.21 0.21 0.53 0.04 −0.36 0.55 −0.04 0.94
PpyrOR59PSE 0.30 0.05 0.51 0.02 −0.21 0.57 0.00 1.00
PpyrOR66 0.35 0.02 0.51 0.01 −0.24 0.36 −0.08 0.64
PpyrOR68 0.59 0.01 0.48 0.02 0.20 0.58 0.09 0.69
PpyrOR61 0.42 0.01 0.38 0.03 0.04 0.96 0.00 1.00

Antennae: M v F

How many genes and ORs DE in male antenna vs female antenna?

antenna_DE <- my.DE.data.df.expanded %>%
  filter(DE_mant_v_fant == "Significant")

mant_v_fant_DE_ORs <- my.DE.OR.data.df.expanded %>%
  filter(DE_mant_v_fant == "Significant") %>%
  arrange(desc(Estimate_mant_v_fant), p_val_adj_mant_v_fant) %>%
  select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
  • 2 genes are significantly differentially expressed between male and female antenna.
  • 1 OR, out of 45 that passed the filter, (2.22)% is DE between male and female antenna.
Which genes and ORs are they?
#for ggplot
my.mant_v_fant.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_fant, y = -log10(p_val_adj_mant_v_fant), color = Direction_mant_v_fant)) +
#  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) + 
  labs(color = "Male v female antenna",
       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated"))

#add circles, label ORco and OR6
my.mant_v_fant.plot.volcano.labels <- my.mant_v_fant.plot.volcano +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_fant, y=-log10(p_val_adj_mant_v_fant), fill = Direction_mant_v_fant), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#F28E2B", "#989ca3", "#4E79A7")) +
  guides(fill = "none") +
  theme(legend.position = "none") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -1.7, vjust = 0.2) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), aes(x = Estimate_mant_v_fant, y = -log10(p_val_adj_mant_v_fant), xend = Estimate_mant_v_fant + 2.5, yend = -log10(p_val_adj_mant_v_fant) + 0.03), color = "grey28") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

my.mant_v_fant.plot.volcano.labels

#for plotly
my.mant_v_fant.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_fant, y = -log10(p_val_adj_mant_v_fant), color = Direction_mant_v_fant, text = name)) +
  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) + 
  labs(color = "Male v female antenna",
#       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
        x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?

my.mant_v_fant.plot.circles <- my.mant_v_fant.plot +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_fant, y=-log10(p_val_adj_mant_v_fant), fill = Direction_mant_v_fant, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#F28E2B", "#989ca3", "#4E79A7")) +
  guides(fill = "none")

ggplotly(my.mant_v_fant.plot.circles)

Circle with bold outline = ORs that passed filtering


DE gene list
mant_v_fant_DE_genes <- my.DE.data.df.expanded %>%
  filter(p_val_adj_mant_v_fant < 0.05) %>%
  arrange(desc(Estimate_mant_v_fant), p_val_adj_mant_v_fant) %>%
  select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)

mant_v_fant_DE_genes.top25 <- mant_v_fant_DE_genes %>%
  slice_head(n = 25)

mant_v_fant_DE_genes.top25 |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(name = "Gene",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
Gene
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p
PREDICTED: Photinus pyralis odorant receptor 94a-like (LOC116160183), transcript variant X1, mRNA 3.48 0.00 0.77 0.03 2.71 0.00 −0.01 0.99
PREDICTED: Photinus pyralis uncharacterized LOC116175393 (LOC116175393), transcript variant X1, mRNA 1.21 0.01 0.00 1.00 1.38 0.05 0.16 0.67

DE OR list
mant_v_fant_DE_ORs <- my.DE.OR.data.df.expanded %>%
  filter(DE_mant_v_fant == "Significant") %>%
  arrange(desc(Estimate_mant_v_fant), p_val_adj_mant_v_fant) %>%
  select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)

mant_v_fant_DE_ORs |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(PpyrOR = "OR",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
OR
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p
PpyrOR6 3.48 0.00 0.77 0.03 2.71 0.00 −0.01 0.99

Leg: M v F

How many genes DE in male leg vs female leg?

leg_DE <- my.DE.data.df.expanded %>%
  filter(DE_mleg_v_fleg == "Significant")

mleg_v_fleg_DE_ORs <- my.DE.OR.data.df.expanded %>%
  filter(DE_mleg_v_fleg == "Significant") %>%
  arrange(desc(Estimate_mleg_v_fleg), p_val_adj_mleg_v_fleg) %>%
  select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
  • 10 genes are significantly differentially expressed between male and female legs
  • 0 ORs, out of 45 that passed the filter, (0)% are DE between male and female legs.
Which genes are they?
#for ggplot
my.mleg_v_fleg.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mleg_v_fleg, y = -log10(p_val_adj_mleg_v_fleg), color = Direction_mleg_v_fleg)) +
#  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) + 
  labs(color = "Male v female leg",
       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated"))

#add circles, label ORco and OR6
my.mleg_v_fleg.plot.volcano.labels <- my.mleg_v_fleg.plot.volcano +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mleg_v_fleg, y=-log10(p_val_adj_mleg_v_fleg), fill = Direction_mleg_v_fleg), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#989ca3", "#4E79A7", "#F28E2B")) +
  guides(fill = "none") +
  theme(legend.position = "none") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -2.2, vjust = 0.2) +
  geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), aes(x = Estimate_mleg_v_fleg, y = -log10(p_val_adj_mleg_v_fleg), xend = Estimate_mleg_v_fleg + 3.3, yend = -log10(p_val_adj_mleg_v_fleg) + 0.03), color = "grey28")+
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

my.mleg_v_fleg.plot.volcano.labels

#for plotly
my.mleg_v_fleg.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mleg_v_fleg, y = -log10(p_val_adj_mleg_v_fleg), color = Direction_mleg_v_fleg, text = name)) +
  geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
  geom_point(size = 1) +
  theme_set(theme_classic() +
              theme(
                axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
                axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
                plot.title = element_text(hjust = 0.5)
            )) +
  coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) + 
  labs(color = "Male v female leg",
#       x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
        x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
                     labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?

my.mleg_v_fleg.plot.circles <- my.mleg_v_fleg.plot +
  geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mleg_v_fleg, y=-log10(p_val_adj_mleg_v_fleg), fill = Direction_mleg_v_fleg, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
  scale_fill_manual(values =  c("#989ca3", "#4E79A7", "#F28E2B")) +
  guides(fill = "none")

ggplotly(my.mleg_v_fleg.plot.circles)

DE gene list
mleg_v_fleg_DE_genes <- my.DE.data.df.expanded %>%
  filter(p_val_adj_mleg_v_fleg < 0.05) %>%
  arrange(desc(Estimate_mleg_v_fleg), p_val_adj_mleg_v_fleg) %>%
  select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)

mleg_v_fleg_DE_genes.top25 <- mleg_v_fleg_DE_genes %>%
  slice_head(n = 25)

mleg_v_fleg_DE_genes.top25 |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(name = "Gene",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
Gene
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p
PREDICTED: Photinus pyralis N-acylneuraminate-9-phosphatase (LOC116164423), mRNA −1.36 0.00 −0.20 0.23 0.25 0.63 1.40 0.03
PREDICTED: Photinus pyralis uncharacterized LOC116178726 (LOC116178726), mRNA −0.19 0.33 0.15 0.50 0.78 0.13 1.12 0.04
PREDICTED: Photinus pyralis trifunctional nucleotide phosphoesterase protein YfkN (LOC116158805), transcript variant X1, mRNA −0.43 0.04 −0.18 0.35 0.76 0.10 1.01 0.04
PREDICTED: Photinus pyralis uncharacterized LOC116174788 (LOC116174788), ncRNA −1.02 0.01 −0.10 0.60 0.06 0.95 0.98 0.04
PREDICTED: Photinus pyralis mpv17-like protein (LOC116176952), mRNA −0.33 0.07 0.00 1.00 0.51 0.20 0.84 0.04
PREDICTED: Photinus pyralis titin homolog (LOC116173520), transcript variant X1, mRNA 0.04 0.75 0.22 0.12 0.49 0.13 0.68 0.04
PREDICTED: Photinus pyralis cysteine dioxygenase type 1 (LOC116159334), mRNA 0.01 0.94 −0.79 0.01 0.03 0.97 −0.78 0.04
PREDICTED: Photinus pyralis uncharacterized LOC116174750 (LOC116174750), transcript variant X1, mRNA 0.02 0.87 −0.05 0.64 −1.42 0.09 −1.49 0.04
PREDICTED: Photinus pyralis phospholipase A1-like (LOC116181243), mRNA 0.00 1.00 −3.09 0.01 −0.36 0.87 −3.45 0.04
PREDICTED: Photinus pyralis cytochrome P450 4c21-like (LOC116176575), mRNA −0.31 0.16 −0.99 0.01 −2.78 0.09 −3.47 0.04

DE OR list
mleg_v_fleg_DE_ORs |>
  gt() |>
  tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
  tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
  tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
  tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
  cols_label(PpyrOR = "OR",
             Estimate_mant_v_mleg = "log2FC", 
             p_val_adj_mant_v_mleg = "p",
             Estimate_fant_v_fleg = "log2FC",
             p_val_adj_fant_v_fleg = "p",
             Estimate_mant_v_fant = "log2FC",
             p_val_adj_mant_v_fant = "p",
             Estimate_mleg_v_fleg = "log2FC",
             p_val_adj_mleg_v_fleg = "p") |>
  fmt_number(decimals = 2) |>
  tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
   tab_style(style = cell_fill(color = "#D4D4D4"),
            locations = cells_body(columns = ends_with("mant_v_fant"))) |> 
  tab_style(style = cell_text(weight = "bold"), 
            locations = list(cells_column_labels(), cells_column_spanners(spanners =  c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))
OR
Male tissue AvL
Female tissue AvL
Antenna MvF
Leg MvF
log2FC p log2FC p log2FC p log2FC p

Are upregulated ORs shared or unique?

x = list(mant_v_mleg_DE_ORs$PpyrOR, fant_v_fleg_DE_ORs$PpyrOR, mant_v_fant_DE_ORs$PpyrOR, mleg_v_fleg_DE_ORs$PpyrOR)

ggVennDiagram(x, category.names = c("Antenna v leg, M", "Antenna v leg, F", "Antenna, M v F", "Leg, M v F")) +
  ggplot2::scale_fill_gradient(low  = "#A0CBE8", high = "#4E79A7") + scale_x_continuous(expand = expansion(mult = .2))

Antenna-upregulated ORs shared between males and females

up_ant_in_both <- intersect(mant_v_mleg_DE_ORs$PpyrOR, fant_v_fleg_DE_ORs$PpyrOR) #shared
sort(up_ant_in_both)
##  [1] "Ppyr/Orco"   "PpyrOR11"    "PpyrOR12"    "PpyrOR16"    "PpyrOR18"   
##  [6] "PpyrOR22"    "PpyrOR24"    "PpyrOR38"    "PpyrOR41"    "PpyrOR53"   
## [11] "PpyrOR56PSE" "PpyrOR58"    "PpyrOR6"     "PpyrOR61"    "PpyrOR62"   
## [16] "PpyrOR63"    "PpyrOR65"    "PpyrOR66"    "PpyrOR68"    "PpyrOR69"   
## [21] "PpyrOR70"    "PpyrOR71"    "PpyrOR75"    "PpyrOR76"    "PpyrOR77"   
## [26] "PpyrOR88"    "PpyrOR92"    "PpyrOR93"    "PpyrOR98"    "PpyrOR99"

Antenna-upregulated ORs unique to males

up_ant_in_males <- setdiff(mant_v_mleg_DE_ORs$PpyrOR, fant_v_fleg_DE_ORs$PpyrOR) #male unique
sort(up_ant_in_males)
## [1] "PpyrOR100" "PpyrOR67"

Antenna-upregulated ORs unique to females

up_ant_in_females <- setdiff(fant_v_fleg_DE_ORs$PpyrOR, mant_v_mleg_DE_ORs$PpyrOR) #female unique
sort(up_ant_in_females)
## [1] "PpyrOR15"    "PpyrOR23"    "PpyrOR28"    "PpyrOR59PSE" "PpyrOR64"   
## [6] "PpyrOR81"

Which ORs that are expressed aren’t DE?

not.in.males <- setdiff(my.DE.OR.data.df.expanded$PpyrOR, mant_v_mleg_DE_ORs$PpyrOR)
not.in.females <- setdiff(my.DE.OR.data.df.expanded$PpyrOR, fant_v_fleg_DE_ORs$PpyrOR)

try2<-my.DE.OR.data.df.expanded[which(my.DE.OR.data.df.expanded$PpyrOR %in% intersect(not.in.males, not.in.females)), ]

try2.long <- try2 %>% pivot_longer(cols = starts_with("SEL"),
                                names_to = "Sample",
                                values_to = "Expression") 

try2.longer <- try2.long %>% mutate(Sex = case_when(Sample == "SEL534ant" ~ "M",
                                                                      Sample == "SEL534BL" ~ "M",
                                                                      Sample == "SEL535ant" ~ "M",
                                                                      Sample == "SEL535BL" ~ "M",
                                                                      Sample == "SEL536ant" ~ "M",
                                                                      Sample == "SEL536BL" ~ "M",
                                                                      Sample == "SEL543ant" ~ "M",
                                                                      Sample == "SEL543BL" ~ "M",
                                                                      Sample == "SEL562ant" ~ "F",
                                                                      Sample == "SEL562BL" ~ "F",
                                                                      Sample == "SEL627ant" ~ "F",
                                                                      Sample == "SEL627BL" ~ "F",
                                                                      Sample == "SEL672ant" ~ "F",
                                                                      Sample == "SEL672BL" ~ "F",
                                                                      TRUE ~ "not yet"))

try2.longer$Sex <- factor(try2.longer$Sex, levels = c("M", "F"))

#add tissue column
try2.longer <- try2.longer %>% mutate(Tissue = case_when(Sample == "SEL534ant" ~ "antenna",
                                                                      Sample == "SEL534BL" ~ "leg",
                                                                      Sample == "SEL535ant" ~ "antenna",
                                                                      Sample == "SEL535BL" ~ "leg",
                                                                      Sample == "SEL536ant" ~ "antenna",
                                                                      Sample == "SEL536BL" ~ "leg",
                                                                      Sample == "SEL543ant" ~ "antenna",
                                                                      Sample == "SEL543BL" ~ "leg",
                                                                      Sample == "SEL562ant" ~ "antenna",
                                                                      Sample == "SEL562BL" ~ "leg",
                                                                      Sample == "SEL627ant" ~ "antenna",
                                                                      Sample == "SEL627BL" ~ "leg",
                                                                      Sample == "SEL672ant" ~ "antenna",
                                                                      Sample == "SEL672BL" ~ "leg",
                                                                      TRUE ~ "not yet"))

#all.DE.ORs.idx.long$Tissue<- factor(all.DE.ORs.idx.long$Tissue, levels = c("antenna", "leg"))

#make sex*tissue column
try2.longer <- try2.longer %>% mutate(Treatment = paste(Sex, Tissue, sep = " "))

try2.longer$Treatment<- factor(try2.longer$Treatment, levels = c("M antenna", "F antenna", "M leg", "F leg"))

ggplot(try2.longer, aes(x = interaction(Sex, Tissue), y = Expression)) +
  geom_boxplot(aes(color = Sex), outlier.shape = NA) + #outlier.size = 0.5
  geom_jitter(size = 0.5, color = "black") +
  facet_wrap2(~ PpyrOR, ncol = 2, scales = "free_y") + #strip = strip
  guides(x = "axis_nested") +
#  scale_fill_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
  scale_color_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
#  ylim(c(8,14.8)) +
  ylab("Normalized counts") +
  xlab("Sample") +
  guides(color = "none")

  • It seems like these genes have higher variance than others in the analysis.

  • Conclusions:
    • Many ORs are upregulated in antennae in both males (32) and females (36)
    • PpyrOR6 is a top candidate for being involved in male-specific behavior (upregulated in antennae of both sexes, and upregulated in male antennae compared to females)
    • PpyrOR100 and PpyrOR67 are also upregulated in male antennae vs. leg, but not females (and not sig. DE between sexes in antennae)
    • PpyrOR15, PpyrOR23, PpyrOR28, PpyrOR59PSE, PpyrOR64, PpyrOR81 are upregulated in female antennae vs leg, but not males (and not sig. DE between sexes in antennae)

Testing model adequacy for DE ORs

# First, index the rows of the dataframe so that each gene has an index#
my.DE.data.df.expanded.idx <- my.DE.data.df.expanded %>%
  mutate(index = 1:13340)

#nrow(my.DE.data.df.expanded.idx)

#then extract the DE ORs
all.DE.ORs.idx <- my.DE.data.df.expanded.idx %>%
  filter(!is.na(PpyrOR)) %>%
  filter(DE_mant_v_mleg == "Significant" | DE_fant_v_fleg == "Significant" | DE_mant_v_fant == "Significant" | DE_mleg_v_fleg == "Significant")
#nrow(all.DE.ORs.idx) #38 total DE

#then plot the fits

for(i in all.DE.ORs.idx$index){
  print(plot(res_reb@fits[[i]], main = my.DE.data.df.expanded.idx$PpyrOR[i]))
}

Plots

ORs are generally upregulated in antennae

Distribution of DE estimates

my.DE.data.df.expanded.est.longer <- my.DE.data.df.expanded %>%
  select(name, Estimate_mant_v_mleg, Estimate_fant_v_fleg, Estimate_mant_v_fant, Estimate_mleg_v_fleg, OR) %>%
  pivot_longer(cols = starts_with("Estimate"), names_to = "Comparison", values_to = "Estimate")

my.DE.data.df.expanded.est.longer$Comparison <- factor(my.DE.data.df.expanded.est.longer$Comparison, 
                                                       levels = c("Estimate_mant_v_mleg", "Estimate_fant_v_fleg", "Estimate_mant_v_fant", "Estimate_mleg_v_fleg"), 
                                                       labels = c("M:AvL", "F:AvL", "A:MvF", "L:MvF"))
                          


ggplot(my.DE.data.df.expanded.est.longer, aes(x=Estimate, fill = OR)) +
  geom_histogram(alpha = 0.5, position = "identity") +
  facet_grid(OR ~ Comparison, scales = "free_y") +
  geom_vline(xintercept = 0, color = "gray", linetype = "dashed") +
  guides(fill = "none")

  • Note that there are no negative values for ORs in Antenna v Leg comparisons – the distribution is skewed right.

1:1 plots of leg vs antenna by sex

male.tissue.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=male.ant.AVG, x=male.leg.AVG, color = Direction_mant_v_mleg)) + 
  geom_point(shape=19, size=1, alpha = 0.75) +
  ggtitle("Male") +
  ylab("Antenna normalized counts") +
  xlab("Leg normalized counts") +
  theme_bw() +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Antenna-biased", "Not significant", "Leg-biased")) +
  geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
  theme_classic() +
  theme(legend.title = element_blank())

male.tissue.1.2.1.labeled <- male.tissue.1.2.1 +
  geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = male.ant.AVG, x = male.leg.AVG), shape = 21, color = "Black") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2)


female.tissue.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=female.ant.AVG, x=female.leg.AVG, color = Direction_fant_v_fleg)) + 
  geom_point(shape=19, size=1, alpha = 0.75) +
  ggtitle("Female") +
  ylab("Antenna normalized counts") +
  xlab("Leg normalized counts") +
  theme_bw() +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Antenna-biased", "Not significant", "Leg-biased")) +
  geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
  theme_classic() +
  theme(legend.title = element_blank())

female.tissue.1.2.1.labeled <- female.tissue.1.2.1 +
  geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = female.ant.AVG, x = female.leg.AVG), shape = 21, color = "Black") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2)

tissue.1.1.plots <- male.tissue.1.2.1.labeled / female.tissue.1.2.1.labeled +
  plot_layout(guides = "collect", axis_titles = "collect") +
  plot_annotation(tag_levels = 'A')

tissue.1.1.plots

#ggsave("2025_03_12_121_plot_tissue.png", plot = tissue.1.1.plots)

PpyrOR6 is the only DE OR (one of two DE genes total) between male and female antennae

1:1 plot of male v female antennae (and legs)

ant.mant_v_fant.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=male.ant.AVG, x=female.ant.AVG, color = Direction_mant_v_fant)) + 
  geom_point(shape=19, size=1, alpha = 0.75) +
  ggtitle("Antenna") +
  ylab("Male normalized counts") +
  xlab("Female normalized counts") +
  theme_bw() +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Male-biased", "Not significant", "Female-biased")) +
  geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
  theme_classic() +
  theme(legend.title = element_blank(), legend.position = "none")

ant.mant_v_fant.1.2.1.labeled <- ant.mant_v_fant.1.2.1 +
  geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = male.ant.AVG, x = female.ant.AVG), shape = 21, color = "Black") +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
  geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2)

leg.mleg_v_fleg.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=male.leg.AVG, x=female.leg.AVG, color = Direction_mleg_v_fleg)) + 
  geom_point(shape=19, size=1, alpha = 0.75) +
  ggtitle("Leg") +
  ylab("Male normalized counts") +
  xlab("Female normalized counts") +
  theme_bw() +
  scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Male-biased", "Not significant", "Female-biased")) +
  geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
  theme_classic() +
  theme(legend.title = element_blank())

leg.mleg_v_fleg.1.2.1.labeled <- leg.mleg_v_fleg.1.2.1 +
  geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = male.leg.AVG, x = female.leg.AVG), shape = 21, color = "Black")

ant.mant_v_fant.1.2.1.labeled / leg.mleg_v_fleg.1.2.1.labeled +
  plot_layout(guides = "collect", axis_titles = "collect") +
  plot_annotation(tag_levels = 'A')

  • The dotted line shows the 1:1 line

Volcano plots

volcanos.plot <- (my.mant_v_mleg.plot.volcano.labels + my.fant_v_fleg.plot.volcano.labels) / (my.mant_v_fant.plot.volcano.labels + my.mleg_v_fleg.plot.volcano.labels) +
   plot_annotation(tag_levels = 'a') +
  plot_layout(guides = "collect")

volcanos.plot

#ggsave("2025_03_12_volcanoes.svg", plot = volcanos.plot, device = "svg", width = 300, height = 300, dpi = 300, units = "mm")
  • Top left, male antennae vs leg
  • Top right, female antennae vs leg
  • Bottom left, male antennae vs female antennae
  • Bottom right, male leg vs female leg

Plots of DE OR + Orco expression w/ p-values and sig bars/letters

Boxplot

all.DE.ORs.idx.long <- all.DE.ORs.idx %>% pivot_longer(cols = starts_with("SEL"),
                                names_to = "Sample",
                                values_to = "Expression")
#length(unique(all.DE.ORs.idx.long$PpyrOR))

#add sex column
all.DE.ORs.idx.long <- all.DE.ORs.idx.long %>% mutate(Sex = case_when(Sample == "SEL534ant" ~ "M",
                                                                      Sample == "SEL534BL" ~ "M",
                                                                      Sample == "SEL535ant" ~ "M",
                                                                      Sample == "SEL535BL" ~ "M",
                                                                      Sample == "SEL536ant" ~ "M",
                                                                      Sample == "SEL536BL" ~ "M",
                                                                      Sample == "SEL543ant" ~ "M",
                                                                      Sample == "SEL543BL" ~ "M",
                                                                      Sample == "SEL562ant" ~ "F",
                                                                      Sample == "SEL562BL" ~ "F",
                                                                      Sample == "SEL627ant" ~ "F",
                                                                      Sample == "SEL627BL" ~ "F",
                                                                      Sample == "SEL672ant" ~ "F",
                                                                      Sample == "SEL672BL" ~ "F",
                                                                      TRUE ~ "not yet"))

all.DE.ORs.idx.long$Sex <- factor(all.DE.ORs.idx.long$Sex, levels = c("M", "F"))

#add tissue column
all.DE.ORs.idx.long <- all.DE.ORs.idx.long %>% mutate(Tissue = case_when(Sample == "SEL534ant" ~ "antenna",
                                                                      Sample == "SEL534BL" ~ "leg",
                                                                      Sample == "SEL535ant" ~ "antenna",
                                                                      Sample == "SEL535BL" ~ "leg",
                                                                      Sample == "SEL536ant" ~ "antenna",
                                                                      Sample == "SEL536BL" ~ "leg",
                                                                      Sample == "SEL543ant" ~ "antenna",
                                                                      Sample == "SEL543BL" ~ "leg",
                                                                      Sample == "SEL562ant" ~ "antenna",
                                                                      Sample == "SEL562BL" ~ "leg",
                                                                      Sample == "SEL627ant" ~ "antenna",
                                                                      Sample == "SEL627BL" ~ "leg",
                                                                      Sample == "SEL672ant" ~ "antenna",
                                                                      Sample == "SEL672BL" ~ "leg",
                                                                      TRUE ~ "not yet"))

all.DE.ORs.idx.long$Tissue<- factor(all.DE.ORs.idx.long$Tissue, levels = c("antenna", "leg"))

#make sex*tissue column
all.DE.ORs.idx.long <- all.DE.ORs.idx.long %>% mutate(Treatment = paste(Sex, Tissue, sep = " "))

all.DE.ORs.idx.long$Treatment<- factor(all.DE.ORs.idx.long$Treatment, levels = c("M antenna", "F antenna", "M leg", "F leg"))

#make boxplots
df <- all.DE.ORs.idx.long

df$PpyrOR <- factor(df$PpyrOR, levels = c("Ppyr/Orco", "PpyrOR6", "PpyrOR11", "PpyrOR12", "PpyrOR15", "PpyrOR16", "PpyrOR18", "PpyrOR22", "PpyrOR23", "PpyrOR24", "PpyrOR28",  "PpyrOR38", "PpyrOR41", "PpyrOR53", "PpyrOR56PSE", "PpyrOR58", "PpyrOR59PSE", "PpyrOR61", "PpyrOR62", "PpyrOR63", "PpyrOR64", "PpyrOR65", "PpyrOR66", "PpyrOR67", "PpyrOR68", "PpyrOR69", "PpyrOR70", "PpyrOR71", "PpyrOR75", "PpyrOR76", "PpyrOR77", "PpyrOR81", "PpyrOR88", "PpyrOR92", "PpyrOR93", "PpyrOR98", "PpyrOR99", "PpyrOR100"))

df$facet <- df$PpyrOR

rectangle_fill_sig_both <- "#bdbdbd"
rectangle_fill_sig_male_only <- "#636363"
rectangle_fill_sig_female_only <- "#f0f0f0"

strip <- strip_themed(background_x = elem_list_rect(fill = c(rectangle_fill_sig_both,  #1, Orco           *
                                                             rectangle_fill_sig_both, #2, Or6             *
                                                             rectangle_fill_sig_both,  #3, Or11           *
                                                             rectangle_fill_sig_both, #4, Or12            *
                                                             rectangle_fill_sig_female_only, #5 , Or15    *
                                                             rectangle_fill_sig_both, #6, OR16            *
                                                             rectangle_fill_sig_both,  #7, OR18           *
                                                             rectangle_fill_sig_both, #8, OR22            *
                                                             rectangle_fill_sig_female_only,  #9, OR23    *
                                                             rectangle_fill_sig_both, #10, OR24           *
                                                             rectangle_fill_sig_female_only, #11, OR28    *
                                                             rectangle_fill_sig_both,#12, OR38            *
                                                             rectangle_fill_sig_both, # OR41              *
                                                             rectangle_fill_sig_both,#14, OR53            *
                                                             rectangle_fill_sig_both, # OR56PSE
                                                             rectangle_fill_sig_both,#16, OR58
                                                             rectangle_fill_sig_female_only, #, OR59PSE 
                                                             rectangle_fill_sig_both,#18, OR61
                                                             rectangle_fill_sig_both, # 62
                                                             rectangle_fill_sig_both,#20, 63 
                                                             rectangle_fill_sig_female_only, #, OR64
                                                             rectangle_fill_sig_both,#22, OR65
                                                             rectangle_fill_sig_both, #, OR66
                                                             rectangle_fill_sig_male_only,#24, OR67
                                                             rectangle_fill_sig_both, #, OR68
                                                             rectangle_fill_sig_both,#26, OR69
                                                             rectangle_fill_sig_both, # 70
                                                             rectangle_fill_sig_both,#28, 71
                                                             rectangle_fill_sig_both, # OR75
                                                             rectangle_fill_sig_both,#30, 76
                                                             rectangle_fill_sig_both, #, OR77
                                                             rectangle_fill_sig_female_only, #32, OR81
                                                             rectangle_fill_sig_both, #, OR88
                                                             rectangle_fill_sig_both,#34, 92
                                                             rectangle_fill_sig_both, #, 93
                                                             rectangle_fill_sig_both, #36, 98
                                                             rectangle_fill_sig_both,  #99
                                                             rectangle_fill_sig_male_only), #38, 100
                                                    color = c(rep("black", 38)), 
                                                    linewidth = c(1, 2, rep(1,36))), 
                      text_x = elem_list_text(color = c(rep("black", 23), "white", rep("black", 13), "white")))
#free scales
anno.free.Orco <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(14.2, 13.4), 
                   y2 = c(14.5, 13.7), 
                   xstar = c(2, 3), 
                   ystar = c(14.6, 13.8), 
                   lab = c("*", "*"),
                   PpyrOR = c("Ppyr/Orco", "Ppyr/Orco"))


anno.free.OR6 <- data.frame(x1 = c(1, 1, 2), #done
                   x2 = c(3, 2, 4), 
                   y1 = c(13, 12.5, 12.0), 
                   y2 = c(13.3, 12.8, 12.3), 
                   xstar = c(2, 1.5, 3), 
                   ystar = c(13.4, 12.9, 12.4), 
                   lab = c("*", "*", "*"),
                   PpyrOR = c("PpyrOR6", "PpyrOR6", "PpyrOR6"))

anno.free.OR11 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(10.15, 9.9), 
                   y2 = c(10.25, 10.0), 
                   xstar = c(2, 3), 
                   ystar = c(10.3, 10.05), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR11", "PpyrOR11"))

anno.free.OR12 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.26, 9.15), 
                   y2 = c(9.3, 9.2), 
                   xstar = c(2, 3), 
                   ystar = c(9.31, 9.22), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR12", "PpyrOR12"))

anno.free.OR15 <- data.frame(x1 = c(2), #done, female only
                   x2 = c(4), 
                   y1 = c(9.6), 
                   y2 = c(9.65), 
                   xstar = c(3), 
                   ystar = c(9.67), 
                   lab = c("*"),
                   PpyrOR = c("PpyrOR15"))

anno.free.OR16 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.72, 9.6), 
                   y2 = c(9.77, 9.65), 
                   xstar = c(2, 3), 
                   ystar = c(9.79, 9.67), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR16", "PpyrOR16"))

anno.free.OR18 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(10.3, 10.1), 
                   y2 = c(10.4, 10.2), 
                   xstar = c(2, 3), 
                   ystar = c(10.45, 10.22), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR18", "PpyrOR18"))

anno.free.OR22 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.3, 9.1), 
                   y2 = c(9.4, 9.2), 
                   xstar = c(2, 3), 
                   ystar = c(9.45, 9.22), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR22", "PpyrOR22"))

anno.free.OR23 <- data.frame(x1 = c(2), #done, female only
                   x2 = c(4), 
                   y1 = c(10.4), 
                   y2 = c(10.5), 
                   xstar = c(3), 
                   ystar = c(10.6), 
                   lab = c("*"),
                   PpyrOR = c("PpyrOR23"))

anno.free.OR24 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(10, 9.75), 
                   y2 = c(10.1, 9.85), 
                   xstar = c(2, 3), 
                   ystar = c(10.15, 9.9), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR24", "PpyrOR24"))

anno.free.OR28 <- data.frame(x1 = c(2), #done, female only
                   x2 = c(4), 
                   y1 = c(10.1), 
                   y2 = c(10.18), 
                   xstar = c(3), 
                   ystar = c(10.2), 
                   lab = c("*"),
                   PpyrOR = c("PpyrOR28"))

anno.free.OR38 <- data.frame(x1 = c(1, 2), 
                   x2 = c(3, 4), 
                   y1 = c(9.3, 9.1), 
                   y2 = c(9.4, 9.2), 
                   xstar = c(2, 3), 
                   ystar = c(9.45, 9.22), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR38", "PpyrOR38"))

anno.free.OR41 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.6, 9.4), 
                   y2 = c(9.7, 9.5), 
                   xstar = c(2, 3), 
                   ystar = c(9.75, 9.53), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR41", "PpyrOR41"))

anno.free.OR53 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.6, 9.4), 
                   y2 = c(9.7, 9.5), 
                   xstar = c(2, 3), 
                   ystar = c(9.75, 9.53), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR53", "PpyrOR53"))

anno.free.OR56 <- data.frame(x1 = c(1, 2), #done, both
                   x2 = c(3, 4), 
                   y1 = c(9.45, 9.2), 
                   y2 = c(9.50, 9.25), 
                   xstar = c(2, 3), 
                   ystar = c(9.53, 9.28), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR56PSE", "PpyrOR56PSE"))

anno.free.OR58 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.4, 9.25), 
                   y2 = c(9.45, 9.3), 
                   xstar = c(2, 3), 
                   ystar = c(9.48, 9.33), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR58", "PpyrOR58"))

anno.free.OR59 <- data.frame(x1 = c(2), #female only 
                   x2 = c(4), 
                   y1 = c(9), 
                   y2 = c(9.05), 
                   xstar = c(3), 
                   ystar = c(9.065), 
                   lab = c("*"),
                   PpyrOR = c("PpyrOR59PSE"))

anno.free.OR61 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9, 8.9), 
                   y2 = c(9.05, 8.95), 
                   xstar = c(2, 3), 
                   ystar = c(9.065, 8.965), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR61", "PpyrOR61"))

anno.free.OR62 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.6, 9.4), 
                   y2 = c(9.7, 9.5), 
                   xstar = c(2, 3), 
                   ystar = c(9.735, 9.535), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR62", "PpyrOR62"))

anno.free.OR63 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.6, 9.45), 
                   y2 = c(9.65, 9.5), 
                   xstar = c(2, 3), 
                   ystar = c(9.7, 9.535), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR63", "PpyrOR63"))

anno.free.OR64 <- data.frame(x1 = c(2), #done, female only
                   x2 = c(4), 
                   y1 = c(9.25), 
                   y2 = c(9.3), 
                   xstar = c(3), 
                   ystar = c(9.335), 
                   lab = c("*"),
                   PpyrOR = c("PpyrOR64"))

anno.free.OR65 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(10.0, 9.8), 
                   y2 = c(10.1, 9.9), 
                   xstar = c(2, 3), 
                   ystar = c(10.135, 9.935), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR65", "PpyrOR65"))

anno.free.OR66 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.22, 9.1), 
                   y2 = c(9.27, 9.15), 
                   xstar = c(2, 3), 
                   ystar = c(9.3, 9.18), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR66", "PpyrOR66"))

anno.free.OR67 <-  data.frame(x1 = c(1), #done, male only
                   x2 = c(3), 
                   y1 = c(9.1), 
                   y2 = c(9.15), 
                   xstar = c(2), 
                   ystar = c(9.165), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR67"))

anno.free.OR68 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.3, 9.1), 
                   y2 = c(9.4, 9.2), 
                   xstar = c(2, 3), 
                   ystar = c(9.45, 9.22), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR68", "PpyrOR68"))

anno.free.OR69 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.8, 9.55), 
                   y2 = c(9.9, 9.65), 
                   xstar = c(2, 3), 
                   ystar = c(9.95, 9.7), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR69", "PpyrOR69"))

anno.free.OR70 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(10.1, 9.8), 
                   y2 = c(10.2, 9.9), 
                   xstar = c(2, 3), 
                   ystar = c(10.25,9.95), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR70", "PpyrOR70"))

anno.free.OR71 <- data.frame(x1 = c(1, 2), 
                   x2 = c(3, 4), 
                   y1 = c(9.45, 9.3), 
                   y2 = c(9.5, 9.35), 
                   xstar = c(2, 3), 
                   ystar = c(9.535, 9.385), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR71", "PpyrOR71"))

anno.free.OR75 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.55, 9.4), 
                   y2 = c(9.6, 9.45), 
                   xstar = c(2, 3), 
                   ystar = c(9.635, 9.485), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR75", "PpyrOR75"))

anno.free.OR76 <- data.frame(x1 = c(1, 2), 
                   x2 = c(3, 4), 
                   y1 = c(9.6, 9.4), 
                   y2 = c(9.65, 9.45), 
                   xstar = c(2, 3), 
                   ystar = c(9.7, 9.5), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR76", "PpyrOR76"))

anno.free.OR77 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.3, 9.1), 
                   y2 = c(9.4, 9.2), 
                   xstar = c(2, 3), 
                   ystar = c(9.45, 9.22), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR77", "PpyrOR77"))

anno.free.OR81 <- data.frame(x1 = c(2), #done, female only
                   x2 = c(4), 
                   y1 = c(9.5), 
                   y2 = c(9.55), 
                   xstar = c(3), 
                   ystar = c(9.58), 
                   lab = c("*"),
                   PpyrOR = c("PpyrOR81"))

anno.free.OR88 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(10.3, 10), 
                   y2 = c(10.4, 10.1), 
                   xstar = c(2, 3), 
                   ystar = c(10.5, 10.2), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR88", "PpyrOR88"))

anno.free.OR92 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(10.6, 10.3), 
                   y2 = c(10.7, 10.4), 
                   xstar = c(2, 3), 
                   ystar = c(10.8, 10.5), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR92", "PpyrOR92"))

anno.free.OR93 <- data.frame(x1 = c(1, 2), #done
                   x2 = c(3, 4), 
                   y1 = c(9.75, 9.6), 
                   y2 = c(9.8, 9.65), 
                   xstar = c(2, 3), 
                   ystar = c(9.83, 9.68), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR93", "PpyrOR93"))

anno.free.OR98 <- data.frame(x1 = c(1, 2), #done 
                   x2 = c(3, 4), 
                   y1 = c(9.75, 9.6), 
                   y2 = c(9.8, 9.65), 
                   xstar = c(2, 3), 
                   ystar = c(9.83, 9.68), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR98", "PpyrOR98"))

anno.free.OR99 <- data.frame(x1 = c(1, 2), 
                   x2 = c(3, 4), 
                   y1 = c(10.1, 9.8), 
                   y2 = c(10.2, 9.9), 
                   xstar = c(2, 3), 
                   ystar = c(10.3, 10.0), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR99", "PpyrOR99"))

anno.free.OR100 <- data.frame(x1 = c(1), #done, male only
                   x2 = c(3), 
                   y1 = c(9.45), 
                   y2 = c(9.5), 
                   xstar = c(2), 
                   ystar = c(9.535), 
                   lab = c("*", "*"),
                   PpyrOR = c("PpyrOR100"))

anno.free <- rbind(anno.free.Orco, anno.free.OR6, anno.free.OR11, anno.free.OR12, anno.free.OR15, anno.free.OR16, anno.free.OR18, anno.free.OR22, anno.free.OR23, anno.free.OR24, anno.free.OR28, anno.free.OR38, anno.free.OR41, anno.free.OR53, anno.free.OR56, anno.free.OR58, anno.free.OR59, anno.free.OR61, anno.free.OR62, anno.free.OR63, anno.free.OR64, anno.free.OR65, anno.free.OR66, anno.free.OR67, anno.free.OR68, anno.free.OR69, anno.free.OR70, anno.free.OR71, anno.free.OR75, anno.free.OR76, anno.free.OR77, anno.free.OR81, anno.free.OR88, anno.free.OR92, anno.free.OR93, anno.free.OR98, anno.free.OR99, anno.free.OR100)

anno.free$PpyrOR <- factor(anno.free$PpyrOR, levels = c("Ppyr/Orco", "PpyrOR6", "PpyrOR11", "PpyrOR12", "PpyrOR15", "PpyrOR16", "PpyrOR18", "PpyrOR22", "PpyrOR23", "PpyrOR24", "PpyrOR28",  "PpyrOR38", "PpyrOR41", "PpyrOR53", "PpyrOR56PSE", "PpyrOR58", "PpyrOR59PSE", "PpyrOR61", "PpyrOR62", "PpyrOR63", "PpyrOR64", "PpyrOR65", "PpyrOR66", "PpyrOR67", "PpyrOR68", "PpyrOR69", "PpyrOR70", "PpyrOR71", "PpyrOR75", "PpyrOR76", "PpyrOR77", "PpyrOR81", "PpyrOR88", "PpyrOR92", "PpyrOR93", "PpyrOR98", "PpyrOR99", "PpyrOR100"))

                   
boxplots.free <-ggplot(df, aes(x = interaction(Sex, Tissue), y = Expression)) +
  geom_boxplot(aes(color = Sex), outlier.shape = NA) + #outlier.size = 0.5
  geom_jitter(size = 0.5, color = "black") +
  facet_wrap2(~ PpyrOR, ncol = 5, strip = strip, scales = "free_y") +
  guides(x = "axis_nested") +
#  scale_fill_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
  scale_color_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
 # ylim(c(8,14.8)) +
  ylab("Normalized counts") +
  xlab("Sample") +
  guides(color = "none") +
  geom_text(data = anno.free, aes(x = xstar,  y = ystar, label = lab)) +
  geom_segment(data = anno.free, aes(x = x1, xend = x1, 
           y = y1, yend = y2),
           colour = "black") +
  geom_segment(data = anno.free, aes(x = x2, xend = x2, 
           y = y1, yend = y2),
           colour = "black") +
  geom_segment(data = anno.free, aes(x = x1, xend = x2, 
           y = y2, yend = y2),
           colour = "black") +
  xlab(element_blank())

boxplots.free

#ggsave("2025_03_12_boxplots_free.svg", plot = boxplots.free, device = "svg", width = 250, height = 350, dpi = 400, units = "mm")
  • * = p < 0.05
  • blue = male, orange = female
  • black jittered points = values for each sample
  • shading of heading means significantly DE in male antennae vs legs (dark gray), female antennae v. legs (light gray), or both (medium gray)
  • bold outline means significantly DE in male vs female antennae

Fold-change on linkage groups

#read in data file with this info
my.LG.data <- read.table("./bigger_OR_tibble_2025_03_12.txt", sep = "\t", header = TRUE)

#extract just one total length entry per OR
only.totals.tib <- my.LG.data %>% filter(Feature_title == "Total length")
#length(only.totals.tib$OR) #101

#make orco name compatible
#length(unique(my.DE.data.df.expanded$PpyrOR)) #46 - the ones that passed filtering
my.DE.data.df.expanded$PpyrOR <- gsub("Ppyr\\/Orco", "Ppyr\\\\Orco", my.DE.data.df.expanded$PpyrOR)

#join with expression dataframe
my.joined.tib <- left_join(only.totals.tib, my.DE.data.df.expanded, by = join_by(OR == PpyrOR))

#factor appropriately
my.joined.tib$OR_GROUP <- factor(my.joined.tib$OR_GROUP, levels = c("Orco", "OR1", "2A", "2B", "3", "4", "6"))

my.joined.tib$LG_common_name <- factor(my.joined.tib$LG_common_name, levels = c("LG1", "LG2", "LG3a (X)", "LG3b", "LG4", "LG5", "LG6", "LG7", "LG8", "LG9", "LG10"))

#my.joined.tib$suffix
my.joined.tib$suffix <-  factor(my.joined.tib$suffix, levels = c("", "NTE", "PSE", "CTE"), labels = c("Full length", "Partial", "PSE", "Partial"))

#plot of LGs
#define option for non-scientific axis notation
options(scipen = 999)

#make tibble of LGs and format it (factor LGs)
LG_tibble <- tibble(LG_common_name = levels(my.joined.tib$LG_common_name), LG_length = c(70905644, 53657282, 22220000, 28664792, 50607672, 49173113, 47017841, 43653355, 37158542, 24427974, 20869423))

LG_tibble$LG_common_name <- factor(LG_tibble$LG_common_name, levels(my.joined.tib$LG_common_name))

my.joined.tib$DE_mant_v_mleg <- factor(my.joined.tib$DE_mant_v_mleg, levels = c("Significant", "Not significant", "NA"))

my.joined.tib$DE_fant_v_fleg <- factor(my.joined.tib$DE_fant_v_fleg, levels = c("Significant", "Not significant", "NA"))

my.joined.tib$DE_mant_v_fant <- factor(my.joined.tib$DE_mant_v_fant, levels = c("Significant", "Not significant", "NA"))

my.joined.tib$DE_mleg_v_fleg <- factor(my.joined.tib$DE_mleg_v_fleg, levels = c("Significant", "Not significant", "NA"))

my.joined.tib$Estimate_mant_v_mleg <- as.numeric(my.joined.tib$Estimate_mant_v_mleg)

my.joined.tib$Estimate_fant_v_fleg <- as.numeric(my.joined.tib$Estimate_fant_v_fleg)

my.joined.tib$Estimate_mant_v_fant <- as.numeric(my.joined.tib$Estimate_mant_v_fant)

my.joined.tib$Estimate_mleg_v_fleg <- as.numeric(my.joined.tib$Estimate_mleg_v_fleg)

library(scales)

#plot LG position
test <- my.joined.tib %>% select(OR, LG_common_name, converted_midpoint, Estimate_mant_v_mleg, Estimate_fant_v_fleg, Estimate_mant_v_fant, Estimate_mleg_v_fleg, DE_mant_v_mleg, DE_fant_v_fleg, DE_mant_v_fant, DE_mleg_v_fleg) %>%
  pivot_longer(
    cols = starts_with("Estimate") | starts_with("DE"),
    cols_vary = "slowest",
    names_to = c(".value", "Comparison"),
    names_pattern = "(.+)_(...._v_....)")

#ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
#  geom_bar(stat = "identity", color = "black", fill = "white") +
#  theme_classic() +
#  coord_flip() +
#  scale_x_discrete(limits = rev) +
#  scale_y_continuous(labels = scales::comma) +
# ylab("Length (bp)") +
#    theme(axis.title.y = element_blank()) +
#  geom_point(data = test, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate, color = DE, size = DE), shape = 21, alpha = 0.7, position = position_jitter(set.seed(42))) +
#  scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
#  scale_color_manual(values = c("black", "NA"), na.value = "gray48") +
#  scale_size_manual(values = c(2, 2), na.value = 0.5) +
#  facet_wrap(~ Comparison)

ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
    theme(axis.title.y = element_blank()) +
    geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_mleg, color = DE_mant_v_mleg, size = DE_mant_v_mleg), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
#  geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_mleg, color = DE_mant_v_mleg, size = DE_mant_v_mleg), width = 0.4, height = 0, shape = 21, alpha = 0.7) +
  scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
  scale_color_manual(values = c("black", "gray48"), na.value = "gray70", name = "", labels = c("adj p \u003c 0.05", "adj p \u2265 0.05", "failed filtering")) +
  scale_size_manual(values = c(2, 1), na.value = 0.5, name = "", labels = c("adj p \u003c 0.05", "adj p \u2265 0.05", "failed filtering")) +
#  theme(legend.position = "none") +
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Mant_v_mleg_fold-change_no_legend_on_LG_2025_03_12.svg", device = "svg")
  
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
    theme(axis.title.y = element_blank()) +
    geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_fant_v_fleg, color = DE_fant_v_fleg, size = DE_fant_v_fleg), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
#  geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_fant_v_fleg), size =2 , width = 0.4, height = 0, shape = 21, alpha = 0.6) +
  scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC Female antennae v legs"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
  theme(legend.position = "none") +
  scale_color_manual(values = c("black", "gray48"), na.value = "gray70") +
  scale_size_manual(values = c(2, 1), na.value = 0.5) +
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Fant_v_fleg_fold-change_on_LG_2025_03_12.svg", device = "svg")

ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
    theme(axis.title.y = element_blank()) +
   geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_fant, color = DE_mant_v_fant, size = DE_mant_v_fant), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
#  geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_fant), size =2 , width = 0.4, height = 0, shape = 21, alpha = 0.6) +
  scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC Male v female antennae"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
  theme(legend.position = "none") +
  scale_color_manual(values = c("black", "gray48"), na.value = "gray70") +
  scale_size_manual(values = c(2, 1), na.value = 0.5) +
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Mant_v_fant_fold-change_on_LG_2025_03_12.svg", device = "svg")

ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
    theme(axis.title.y = element_blank()) +
  geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mleg_v_fleg, color = DE_mleg_v_fleg, size = DE_mleg_v_fleg), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
#geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mleg_v_fleg), size =2 , width = 0.4, height = 0, shape = 21, alpha = 0.6) +
  scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC Male v female legs"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
  theme(legend.position = "none") +
  scale_color_manual(values = c("black", "gray48"), na.value = "gray70") +
  scale_size_manual(values = c(2, 1), na.value = 0.5) +
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Mleg_v_fleg_fold-change_on_LG_2025_03_12.svg", device = "svg")
#for some reason, jittering changes when you assign the plot to a variable (like the set.seed doesn't work - so saved individually and combined in adobe)
  • Top to bottom:
    • Male antennae v leg
    • Female antennae v leg
    • Male antennae v female antennae
    • Male leg v female leg

Expression on linkage groups

#plot of LGs
#define option for non-scientific axis notation
options(scipen = 999)

#extract data for expression
expression_on_LGs.df <- my.joined.tib %>% select(OR, LG_common_name, converted_midpoint, male.ant.AVG, female.ant.AVG, male.leg.AVG, female.leg.AVG, DE_mant_v_mleg, DE_fant_v_fleg, DE_mant_v_fant, DE_mleg_v_fleg) %>%
  pivot_longer(
    cols = ends_with("AVG"),
    cols_vary = "slowest", names_to = "Tissue", values_to = "Expression")

expression_on_LGs.df$Tissue <- factor(expression_on_LGs.df$Tissue, levels = c("male.ant.AVG", "female.ant.AVG", "male.leg.AVG", "female.leg.AVG"), labels = c("Male antennae", "Female antennae", "Male leg", "Female leg"))

expression_on_LGs.df$filtered <- "PASS"
expression_on_LGs.df$filtered[which(is.na(expression_on_LGs.df$Expression))] <- "FAIL"
expression_on_LGs.df$filtered <- factor(expression_on_LGs.df$filtered, levels = c("PASS", "FAIL"), labels = c("", "Filtered out"))

library(scales)
#plot LG position

ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
  theme(axis.title.y = element_blank()) +
  geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Male antennae") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
  scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
  scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
    force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Male_ant_norm_counts_on_LG_2025_03_12.svg", device = "svg")
  
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
  theme(axis.title.y = element_blank()) +
  geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Female antennae") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
  scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
  scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
  theme(legend.position = "none") +
  force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Female_ant_norm_counts_on_LG_2025_03_12.svg", device = "svg")

ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
  theme(axis.title.y = element_blank()) +
  geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Male leg") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
  scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
  scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
  theme(legend.position = "none") +
  force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Male_leg_norm_counts_on_LG_2025_03_12.svg", device = "svg")

ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
  geom_bar(stat = "identity", color = "black", fill = "white") +
  theme_classic() +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  scale_y_continuous(labels = scales::comma) +
  ylab("Length (bp)") +
  theme(axis.title.y = element_blank()) +
  geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Female leg") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
  scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
  scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
  theme(legend.position = "none") +
  force_panelsizes(rows = unit(85, "mm"),
                   cols = unit(85, "mm"))

#ggsave("Female_leg_norm_counts_on_LG_2025_03_12.svg", device = "svg")
#for some reason, jittering changes when you assign the plot to a variable (like the set.seed doesn't work - so saved individually and combined in adobe)
  • Top to bottom:
    • Male antennae
    • Female antennae
    • Male leg
    • Female leg

Expression on phylogeny

#load tree
#read in tree
tree <- read.tree("./run_2.contree_rooted_2025_03_11.newick")

#trim tree
#get the data from the tree
t_data <- as_tibble(tree)

#rename PpyrOR8
t_data <- t_data %>% mutate(label2 = label)
t_data$label2[which(t_data$label == "PpyrOR8b")] <- "PpyrOR8"

#rename Orco
t_data$label2[which(t_data$label == "Ppyr_Orco")] <- "PpyrOrco"

tree.rename <- rename_taxa(tree, t_data, label, label2)
t.rename_data <- as_tibble(tree.rename)

#just get the total OR info, 1 per OR, only firefly
only_totals <- my.joined.tib %>% 
  filter(feature == "total") %>%
  filter(SPECIES == "Photinus pyralis")

#join the data
t_data <- left_join(t.rename_data, only_totals, by = join_by(label == OR_NAME))

#drop tips that are not Ppyr
to_drop <- which(is.na(t_data$LG))
tree_reduced <- drop.tip(tree.rename, to_drop)

#check tree looks good
#ggtree(tree_reduced) + geom_tiplab(size = 2) #good

#get reduced tree data
tr_data <- as_tibble(tree_reduced)

#add in OR info
tr_data <- left_join(tr_data, only_totals, by = join_by(label == OR_NAME))

tr_data$OR <- factor(tr_data$OR, levels = tr_data$OR)

tr_data$OR_GROUP <- factor(tr_data$OR_GROUP, levels = c("Orco", "OR1", "2A", "2B", "3", "4", "6")) 

tr_data$LG_common_name <- factor(tr_data$LG_common_name, levels = c("LG1", "LG2", "LG3a (X)", "LG3b", "LG4", "LG5", "LG6", "LG7", "LG8", "LG9", "LG10"))

expression.tibble <- tr_data %>% select(parent, node, branch.length, label, male.ant.AVG, female.ant.AVG, male.leg.AVG, female.leg.AVG)

expression.tibble.longer <- expression.tibble %>% pivot_longer(ends_with("AVG"), names_to = "Sex x Tissue", values_to = "Normalized counts")

expression.tibble.longer <- expression.tibble.longer %>% mutate(Sex = ifelse(str_detect(`Sex x Tissue`, "female"), "Female", "Male"),
                                                                             Tissue = ifelse(str_detect(`Sex x Tissue`, "ant"), "Antennae", "Legs"))
  
expression.tibble.longer$Sex <- factor(expression.tibble.longer$Sex, levels = c("Male", "Female"), labels = c("M", "F"))

expression.tibble.longer$Tissue <- factor(expression.tibble.longer$Tissue, levels = c("Antennae", "Legs"))

#orient tree
#attach expression data
#plot - can I have a series of circles @ end with expression values instead of 4 different trees? Or a heatmap?
library(aplot)

p <- ggtree(tree_reduced) %<+% tr_data

#add tip labels, group designations (color), and node labels to see which nodes to flip
#p + geom_tiplab(size = 2) + 
#  geom_tippoint(aes(fill = OR_GROUP), color = "black", shape = 21, size = 2) +
#  geom_text(aes(label = node), size = 2, vjust = -1, hjust = 1.2)
#  scale_fill_manual(values = c("#4E79A7","#E15759", "#B6992D", "#F28E2B", "#B07AA1", "#59A14F", "#797067"))

#flip the 2B node to be next to 2A for ease of visualization
p2 <- flip(p, 124, 107)

#plot flipped tree with labels and colors to check 
#p2 + 
#  geom_tiplab(size = 2, hjust = -0.1) + 
#  geom_text(aes(label = node), size = 2, vjust = -1, hjust = 1.2) +
#  geom_tippoint(aes(color = OR_GROUP),  size = 1.5, alpha = 0.9)# +
#  scale_color_manual(values = c("#4E79A7","#E15759", "#B6992D", "#F28E2B", "#B07AA1", "#59A14F", "#797067")) +
#  theme(legend.position = "right") +
# xlim(0,10) +
#  ylim(0, 101)

#get node values for labels
Group2A.MRCA <- getMRCA(tree_reduced, c("PpyrOR6", "PpyrOR2")) #group2A
Group2B.MRCA <- getMRCA(tree_reduced, c("PpyrOR81", "PpyrOR65")) #group2B
Group3.MRCA <-getMRCA(tree_reduced, c("PpyrOR92", "PpyrOR90")) #group3
Group4.MRCA <- getMRCA(tree_reduced, c("PpyrOR84", "PpyrOR87")) #group4
Group6.MRCA <- getMRCA(tree_reduced, c("PpyrOR100", "PpyrOR99")) #group6

z<-p2 +
  geom_tiplab(align=TRUE, color = "white") + hexpand(0.15) + xlim(0,10) + 
  geom_cladelab(node=Group2A.MRCA, label="2A", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
geom_cladelab(node=Group2B.MRCA, label="2B", align = TRUE, offset=0.9, barsize=1, fontsize=7) + 
  geom_cladelab(node=Group3.MRCA, label="3", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
  geom_cladelab(node=Group4.MRCA, label="4", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
  geom_cladelab(node=Group6.MRCA, label="6", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
  geom_strip("PpyrOR1", "PpyrOR1", label="OR1", align = TRUE, offset=-1.0, barsize=1, fontsize=3) +
  annotate("segment", x = 3.45, y = 8.7, xend = 4.2, yend = 8.7, color = "gray48", lty = "1111", linewidth = 0.5)

z2 <- ggplot(expression.tibble.longer, aes(x = interaction(Sex, Tissue), y = label)) +
  geom_tile(aes(fill = `Normalized counts`), na.rm = TRUE) + 
  scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) + 
  theme_classic() + 
  xlab(NULL) + 
  ylab(NULL) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  guides(x = "axis_nested")
#+
#  force_panelsizes(rows = unit(85, "mm"),
#                   cols = unit(45, "mm"))



z2 %>% insert_left(z)

#ggsave("Expression_on_phylo_2025_03_12.svg", device = "svg")

Session information

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.7.4
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] aplot_0.2.3                 scales_1.3.0               
##  [3] stringr_1.5.1               treeio_1.28.0              
##  [5] ape_5.8                     ggtree_3.12.0              
##  [7] tidyr_1.3.1                 plotly_4.10.4              
##  [9] dplyr_1.1.4                 REBEL_0.1.0                
## [11] ggpattern_1.1.1             patchwork_1.3.0.9000       
## [13] ggh4x_0.2.8                 ggplot2_3.5.1              
## [15] ggVennDiagram_1.5.2         gt_0.11.1                  
## [17] readr_2.1.5                 SummarizedExperiment_1.34.0
## [19] Biobase_2.64.0              GenomicRanges_1.56.2       
## [21] GenomeInfoDb_1.40.1         IRanges_2.38.1             
## [23] S4Vectors_0.42.1            BiocGenerics_0.50.0        
## [25] MatrixGenerics_1.16.0       matrixStats_1.4.1          
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        viridisLite_0.4.2       farver_2.1.2           
##  [4] fastmap_1.2.0           lazyeval_0.2.2          digest_0.6.37          
##  [7] lifecycle_1.0.4         tidytree_0.4.6          magrittr_2.0.3         
## [10] compiler_4.4.1          rlang_1.1.4             sass_0.4.9             
## [13] tools_4.4.1             yaml_2.3.10             data.table_1.16.4      
## [16] knitr_1.49              labeling_0.4.3          S4Arrays_1.4.1         
## [19] htmlwidgets_1.6.4       bit_4.5.0.1             DelayedArray_0.30.1    
## [22] xml2_1.3.6              abind_1.4-8             withr_3.0.2            
## [25] purrr_1.0.2             grid_4.4.1              colorspace_2.1-1       
## [28] cli_3.6.3               rmarkdown_2.29          crayon_1.5.3           
## [31] generics_0.1.3          rstudioapi_0.17.1       httr_1.4.7             
## [34] tzdb_0.4.0              cachem_1.1.0            zlibbioc_1.50.0        
## [37] parallel_4.4.1          ggplotify_0.1.2         XVector_0.44.0         
## [40] vctrs_0.6.5             yulab.utils_0.1.8       Matrix_1.7-0           
## [43] jsonlite_1.8.9          gridGraphics_0.5-1      hms_1.1.3              
## [46] bit64_4.5.2             crosstalk_1.2.1         jquerylib_0.1.4        
## [49] glue_1.8.0              stringi_1.8.4           gtable_0.3.5           
## [52] UCSC.utils_1.0.0        munsell_0.5.1           tibble_3.2.1           
## [55] pillar_1.10.1           htmltools_0.5.8.1       GenomeInfoDbData_1.2.12
## [58] R6_2.5.1                vroom_1.6.5             evaluate_1.0.3         
## [61] lattice_0.22-6          ggfun_0.1.8             bslib_0.8.0            
## [64] Rcpp_1.0.14             SparseArray_1.4.8       nlme_3.1-166           
## [67] xfun_0.50               fs_1.6.5                pkgconfig_2.0.3